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In  this  thesis  we  study  the  turbulent  mixing  and  turbulent  combustion 
in  a  model  scramjet  combustor  with  a  Large  Eddy  Simulation  (LES)  strategy. 
LES  resolves  the  large  and  energetic  motions  while  the  small  subscale  motions 
are  modeled.  Here  the  filtered  Navier-Stokes  equations  are  solved  by  a  fifth- 
order  finite  difference  Weighted- Essentially  Non-Oscillatory  (WENO)  scheme 
dimension  by  dimension.  Subgrid  terms  are  closed  by  the  dynamic  Smagorin- 
sky  model.  Chemical  source  terms  are  calculated  directly  using  a  finite  rate 
chemistry  model  with  a  reduced  chemistry  mechanism.  The  equilibrium  tur¬ 
bulent  boundary  layer  model  of  J.  Larsson  is  used  to  calculate  the  shear  stress 


and  heat  flux  at  the  wall.  Inflow  turbulent  is  generated  by  the  digital  filtering 
method. 

The  main  result  is  a  methodology  to  predict  the  mesh  convergence  for 
three-dimensional  turbulent  combustion  simulation,  based  on  a  less  expensive 
suites  of  one- dimensional  and  two-dimensional  simulations.  We  first  deter¬ 
mine  the  grid  requirements  for  finite  rate  chemistry  with  detailed  and  reduced 
chemical  mechanism  respectively  in  the  context  of  one-dimensional  simula¬ 
tions.  These  criteria  are  verified  through  simulation  in  a  two-dimensional 
context  and  refined  with  corrections  due  to  turbulent  transport.  They  are 
then  applied  to  three-dimensional  simulations.  A  grid  sensitivity  study  of  the 
turbulent  boundary  layer  is  conducted  in  a  2D  context. 

Simulation  results  are  validated  through  comparison  with  a  simulation  of 
the  same  problem  conducted  by  J.  Larsson,  using  a  different  methodology  and 
by  comparison  to  experiments  performed  at  Stanford  University. 
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Chapter  1 

Introduction 

1.1  Backgound 

A  jet  engine  is  an  air  breathing  engine  which  carries  fuel  on  board  and 
uses  oxygen  ingested  from  the  atmosphere  for  combustion.  There  are  three 
kinds  of  jet  engine:  turbo-machinery  based  jet  engines,  ramjets,  and  supersonic 
combustion  ramjets  (scramjets).  The  specific  impulses  of  these  jet  engines  are 
shown  in  Fig  1.1.  This  figure  indicates  the  advantage  of  scramjet  when  the 
Mach  number  is  above  six.  The  scramjet  is  the  most  promising  engine  for 
hypersonic  flight.  It  has  been  a  very  hot  research  area  for  last  several  decades. 

The  schematic  plot  of  HyShot-II,  a  scramjet  combustion  chamber,  is 
shown  in  Fig  1.2.  The  structure  of  the  scramjet  combustion  chamber  is  simple. 
The  front  part  (inlet  ramp)  of  the  engine  generates  the  shocks  which  compress 
the  inflow  supersonic  air.  The  velocity  of  the  inflow  air  decreases  as  a  result 
of  the  compression  but  it  is  still  supersonic.  Meanwhile,  the  fuel  is  injected 
(transversely  in  HyShot  II)  into  this  high  enthalpy  crossflow.  After  the  fuel 
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Figure  1.1:  Specific  impulses  of  three  kinds  of  jets  and  rocket  [60]. 


Figure  1.2:  A  schematic  diagram  for  HyShot-II  Model  (Image  source  Gardner 
[23]  ). 
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mixes  with  the  air,  it  burns  and  expands  to  the  outlet  and  generates  thrust 
for  the  aircraft.  There  is  no  ignition  device  needed.  Although  the  structure  is 
simple,  the  physical  phenomena  in  the  combustion  chamber  are  complex  due 
to  shock  waves,  chemical  reactions,  turbulent  mixing  and  turbulent  boundary 
layers.  These  effects  are  nonlinear  and  couples  with  each  other.  This  fact 
makes  the  prediction  and  simulation  of  the  scramjet  combustion  difficult. 

The  research  for  scramjet  engine  has  stretched  over  several  decades.  The 
evolution  of  scramjet  development  can  be  found  in  a  review  paper  of  Fry  [17]. 
Many  flight  and  ground  experiments  have  been  conducted.  Among  them,  the 
HyShot  II  scramjet  engine  designed  by  the  University  of  Queensland,  Aus¬ 
tralia,  was  successfully  launched  in  2002.  The  HEG  facility  of  the  German 
Aerospace  Agency  DLR  has  carried  out  HyShot  II  scramjet  ground  tests  [23]. 

This  thesis  focuses  on  the  numerical  simulation  of  a  model  scramjet  com¬ 
bustor  designed  by  Gamba  [21].  Gamba  carried  out  experiments  to  investigate 
the  mixing,  ignition,  and  combustion  in  supersonic  combustion  in  the  compact, 
optically  accessible  scramjet  combustor  model  at  Stanford  6”  Expansion  Tube 
Facility.  Fig.  1.3  shows  the  configuration  plot  of  the  model  scramjet  combustor. 
The  diagnostic  techniques  of  this  experiment  include  the  PLIF  (Planar  Laser 
Induced  Fluorescence)  imaging  of  OH  concentration  on  orthogonal  planes,  the 
Chemiluminescence  imaging  of  OH*  and  the  pressure  measurement  at  the  top 
wall. 

Besides  flight  and  ground  experiments,  numerous  Computational  Fluid 
Dynamic  (CFD)  techniques  have  been  used  to  understand  the  physics  inside 
scramjet  combustion  chamber.  The  fuel  choices  in  these  simulations  include 
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Figure  1.3:  A  diagram  for  the  model  scramjet  combustor  (Image  source  Garnba 
[21])- 


Figure  1.4:  A  schematic  plot  of  the  transverse  jet  in  supersonic  cross  flow 
(Image  source  Garnba  [20]  ). 


hydrocarbon  fuel  and  hydrogen  fuel.  For  simulation  approach,  both  Reynolds 
Average  Navier-Stokes  (RANS)  approach  [14,  18,  44,  49]  and  Large  Eddy  Sim¬ 
ulation  (LES)  approach  [6,  36,  29]  are  applied.  Various  combustion  models  are 
used  with  either  detailed  or  reduced  chemistry  mechanism.  The  most  popular 
method  among  them  is  the  flamelet  method. 

One  classical  flow  structure  in  scramjet  simulation  is  the  transverse  jet  in 
supersonic  cross  flow  (JICF).  Fig.  1.4  is  the  schematic  plot  of  JICF  which  has 
many  flow  features  [20]:  the  vortical  structures  generated  by  the  interaction 
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of  the  jet  with  the  crossflow,  a  bow  shock  produced  by  the  blockage  effect  of 
the  jet,  the  upstream  recirculation  region  caused  by  the  bow-shock-induced 
boundary  layer  separation,  and  another  recirculation  region  downstream  of 
the  injector.  One  important  parameter,  which  determine  the  characteristics 
of  JICF,  is  the  jet-to-crossflow  momentum  flux  ratio  J  defined  as 

j  _  P}etu]et 
Po 


The  jet  exit  Reynolds  number  defined  as 


Pjeth^jet-^/hjet  (1-2) 

is  another  important  parameter.  Here  D  is  the  diameter  of  H2  injector.  Garnba 
[20]  used  the  OH  PL1F  imaging  to  investigate  the  reaction  in  JICF  at  J  = 
5.0.  He  measured  the  thickness  of  the  diffusion  flame  (OH  layer)  near  the 
injector  and  found  the  thickness  of  the  flame  is  between  300  -  900  microns. 
The  Kolmogorov  length  scale  in  the  scramjet  combustor  is  approximately  10 
microns  (the  estimation  of  the  Kolmogorov  length  scale  is  in  Sec.  5.1).  Thus 
the  Kolmogorov  scale  is  much  smaller  than  the  chemical  scale.  Thus  We  have 
the  numerical  capacity  to  resolve  chemistry  while  we  do  not  have  the  capacity 
to  resolve  turbulence.  The  fact  leads  us  to  choose  the  “finite  rate  chemistry” 
method.  We  will  further  discuss  it  in  Chap.  4. 
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1.2  Thesis  Structure 


This  thesis  presents  our  work  of  the  numerical  simulation  of  the  model 
scramjet  combustor  by  Gamba.  ft  is  organized  as  follows. 

The  first  chapter  explains  the  background  of  scramjet  simulation. 

The  second  chapter  includes  the  major  theories,  governing  equations  and 
numerical  methods  used  in  scramjet  simulation.  Basic  knowledge  of  turbulence 
and  three  turbulence  modelling  methodologies  are  introduced.  The  governing 
Naiver-Stokes  equations  of  compressible  fluid  are  covered.  We  explain  the  main 
modules  in  the  compressible  Frontier  code  we  use  to  simulate  model  scramjet 
combustor: 

1.  The  hyperbolic  module:  The  numerical  discretization  scheme  for  solving 
the  hydrodynamics  dimension  by  dimension. 

2.  The  EOS  module:  The  Equation  of  State  and  transport  properties. 

3.  The  parabolic  module:  Molecular  transport  effects  and  subgrid  scaled 
terms. 

4.  Turbulent  inflow  generation. 

The  third  chapter  focuses  on  the  implementation  and  grid  sensitivity 
study  of  the  turbulent  boundary  layer  through  suits  of  two  dimensional  simu¬ 
lation. 

The  fourth  chapter  studies  the  chemical  combustion  in  the  combustion 
chamber.  The  finite  rate  chemistry  mechanism  with  both  detailed  and  re¬ 
duced  chemical  mechanism  is  explored.  Then  grid  sensitivity  study  of  one- 
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dimensional  laminar  flames  and  two-dimensional  diffusional  flames  are  con¬ 


ducted.  The  thickened  flame  model  is  also  explored. 

The  fifth  chapter  shows  the  result  and  analysis  of  scramjet  3D  simula¬ 
tion,  which  applies  the  resolution  predicted  from  the  third  and  fourth  chapter. 
The  scramjet  3D  simulation  result  is  compared  with  the  experiment  and  the 
Larsson’s  simulation. 

The  sixth  chapter  concludes  the  results  in  this  thesis  and  makes  sugges¬ 
tion  for  future  study. 
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Chapter  2 


Turbulence  Theory,  Governing  Equations  and 
Numerical  Methods 

2.1  Introduction  to  Turbulence 

Turbulence  is  “unsteady,  periodic  motion  in  which  all  three  velocity  com¬ 
ponents  fluctuate,  mixing  matter,  momentum,  and  energy.”  Many  flows  in 
nature  are  turbulent,  see  Fig.  2.1.  Turbulent  flows  occur  at  high  Reynolds 
numbers,  when  different  scales  of  eddies  arise  from  the  complex  interaction 
between  the  viscous  terms  and  the  inertia  terms  in  the  momentum  equations 
[19].  Here,  Reynolds  number  is  defined  as  the  ratio  of  inertial  stress  and  vis¬ 
cous  stress  and  can  often  be  expressed  as  uL/u  (u  is  velocity,  L  is  character 
length  and  v  is  the  viscosity  coefficient). 

At  low  Reynolds  numbers,  flows  are  laminar.  Transitions  from  laminar 
flow  to  turbulence  flow  occur  between  Re  =  2000  to  Re  =  13000.  There  are 
several  characteristics  of  turbulence  [64], 


1.  Randomness.  Turbulent  flows  are  chaotic  and  can  not  be  described  by 


(c)  (d) 


Figure  2.1:  Turbulence  in  nature,  (a)  wake  behind  the  plane,  (b)  cloud  un¬ 
der  Kelvin- Helmholtz  instability,  (c)  “The  Great  Wave  Off  Kanagawa”,  (d) 
turbulent  flame. 
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deterministic  methods.  Instead,  turbulent  flows  are  usually  described 
through  statistical  approaches. 

2.  Diffusivity.  The  diffusivity  of  turbulence  speeds  up  the  mixing  ,  the 
momentum,  heat  and  mass  transfer  process.  The  diffusivity  of  turbu¬ 
lent  is  very  important  in  many  applications,  for  example,  the  turbulence 
combustion  process  of  scramjet  chamber. 

3.  Dissipation.  The  viscous  shear  stress  transfers  energy  from  kinetic  energy 
to  the  internal  energy  of  the  flow.  Thus  a  sustained  turbulence  flow 
requires  a  continuous  supply  of  energy. 

The  existence  of  turbulence  has  a  large  influence  on  the  flow.  For  example, 
the  near  wall  turbulence  increases  the  shear  stress  at  the  wall;  the  turbulence 
in  combustion  chamber  speeds  up  the  mixing  and  the  chemical  reactions. 

Although  turbulence  is  an  important  phenomenon,  understanding  turbu¬ 
lence  is  difficult.  Actually,  turbulence  was  considered  by  Nobel  Laureate  and 
Richard  Feynman  to  be  “one  of  the  biggest  outstanding  problems  in  classical 
physics.”  Two  important  theories  of  turbulence  are  the  energy  cascade  and 
the  Kolmogorov  hypotheses. 

The  turbulent  flow  field  contains  eddies  with  different  sizes.  The  large 
eddies  have  most  of  kinetic  energy.  They  are  unstable  and  transfer  the  energy 
to  smaller  eddies.  The  smaller  eddies  follow  similar  processes,  and  transfer 
the  energy  to  even  smaller  eddies.  The  energy  cascade  is  the  process  that  the 
kinetic  energy  that  enters  the  turbulence  at  the  large  scale  eddies  is  transferred 
to  smaller  and  smaller  scale  eddies,  until  the  Reynolds  number  of  these  eddies 
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is  small  enough  and  these  eddies  are  stable.  The  kinetic  energy  of  the  smallest 
scale  eddies  is  dissipated  away  through  viscosity  and  convert  into  heat.  This 
process  is  summarized  by  L.F.  Richardson  [58]  as: 

Big  whirls  have  little  whirls 
Which  feed  on  their  velocity; 

And  little  whirls  have  lesser  whirls, 

And  so  on  to  viscosity 
in  the  molecular  sense. 

We  denote  the  energy  dissipation  rate  at  which  turbulence  kinetic  energy  is 
converted  into  thermal  internal  energy  as  e  = 

The  Kolmogorov’s  hypothesis  is  based  on  three  hypothesis  together  with 
dimensional  arguments.  The  first  hypothesis  of  similarity  states  that  “in  every 
turbulent  flow  at  sufficiently  high  Reynolds  number,  the  statistics  of  the  small 
scale  motions  (/  <  Iei )  have  a  universal  form  that  is  uniquely  determined  by 
e  and  za”  Here  Iei  is  the  length  scale  which  separates  the  large,  energetic, 
geometric  dependent  scales  from  small  universal  scales.  The  size  range  l  <  Iei 
is  called  the  universal  equilibrium  range.  The  range  where  l  >  Iei  is  called  the 
energy  containing  range. 

The  second  hypothesis  of  similarity  states  that  “  in  every  turbulent  flow 
at  sufficiently  high  Reynolds  number,  the  statistics  of  the  motions  of  scale  l 
in  the  range  (0  >  I  »  r)  have  a  universal  form  that  is  uniquely  determined 
by  e  and  independent  of  za”  This  hypothesis  splits  the  universal  equilibrium 
range  into  two  subranges:  the  inertial  subrange  where  motions  are  uniquely 
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Figure  2.2:  Turbulence  length  scales  [57] 


determined  by  e  and  the  dissipation  range  where  motions  are  affected  by  vis- 
cousity.  Different  length  scales  are  related  to  these  three  ranges  of  turbulence, 
summarized  as  Fig  2.2. 

The  large  and  energetic  eddies  have  a  length  scale  of  the  same  order 
with  the  problem  geometry,  denoted  by  l  and  velocity  length  scale  of  uc.  The 
Kolmogorov  length  scale,  velocity  scale  and  time  scale  are: 


/„3\  !/4 

length  scale:  rjk  =  ^ — J 

(2.1) 

velocity  scale:  uv  =  (V3e) 1 ^ 

(2.2) 

/UX1/2 

time  scale:  tv  =  J 

(2.3) 

The  ratios  of  the  Kolmogorov  scales  to  the  integral  scales  are  related  to 
the  Reynolds  number  by 


Re~3/4 

(2.4) 

.Re- 1/4 

(2.5) 

uc 
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rs_/ 


Re-1'2 


(2.6) 


The  intermediate  scales  in  the  inertial  range  of  turbulence  are  called  as  the 
Taylor  scales  A t-  The  ratio  of  Taylor  scales  to  the  Kolmogorov  scales  and  the 
energetic  large  scales  are: 

y  ~  Re~l/2  (2.7) 

—  ~  Rel/i  (2.8) 

V 

The  distribution  of  the  turbulent  kinetic  energy  over  these  three  scales 
can  be  observed  by  a  spectral  analysis.  In  the  spectrum  analysis  we  represent 
different  length  scales  by  their  wave  number  (k)  and  represent  the  velocity  field 
in  Fourier  space.  By  multiplying  the  velocity  at  a  particular  wave  number  with 
its  complex  conjugate  we  get  the  spectrum  analysis  presenting  the  distribution 
of  kinetic  energy  across  ranges  of  wave  numbers. 

Kolmogorov’s  5/3  spectrum  law  states  that  the  energy  spectrum  E(k)  in 
the  inertial  range  follows  E(k)  =  Ce2^3K~5^3  where  e  is  the  energy  dissipation 
rate,  k  is  XXX  and  C  is  the  universal  Kolmogorov  constant.  The  typical 
energy  spectrum  of  all  three  ranges  are  shown  in  Fig.  2.3. 


2.2  Conservation  Laws  and  Governing  Equations 


2.2.1  Navier-Stokes  Equation 


The  governing  equations  for  our  simulation  are  the  reactive  compressible 
Navier-Stokes  equations  for  an  ideal  reactive  gas  [51,  56],  which  describe  the 
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Figure  2.3:  Schematic  plot  of  a  typical  turbulent  energy  spectrum  [47]. 


conservation  of  mass,  momentum  and  energy: 


dp  d 
dt  dxj 


0, 


(2.9) 


d 


dt 


d 


7h '  pUi  1  +  d^ 1  pUiUj  '  +  ^  = 


dp  dp 


v 


dxi 


dxj  ’ 


dE  d[(E  +  p)uj]  d 
dt  dx;  dxi 


Ti’Ui  '  "  dx 


dqi  +  H  . 


(2.10) 


(2,11) 


together  with  the  balance  equation  for  the  mass  fraction  of  each  species  in  the 


mixture: 

+Ai-  (212) 

Here  p  is  the  density,  u  is  the  velocity,  E  is  the  specific  energy,  is  the  mass 
fraction  of  species  i,  and  rhi  is  the  production  rate  of  species  i.  H  is  the  rate  of 
heat  release  from  the  chemical  reactions.  The  viscous  stress  tensor  Pj  derived 


14 


from  Newtonian  fluid  reflects  the  proportionality  between  the  shear  stress  and 
the  rate  of  deformation: 


Tij  =  2/^*- 


+ 


duk 

dxk 


(2.13) 


where  /i  is  the  dynamic  viscosity  coefficient  and  <5  is  the  Kronecker  delta.  The 
heat  flux  q.j  is  calculated  by  Fourier’s  law  of  conduction,  an  empirical-base 
derivation  relates  the  heat  transfer  with  temperature  gradient, 


*  =  - k s;  •  (2'14) 

where  T  is  the  temperature  and  k  is  the  thermal  conductivity  coefficient.  Cal¬ 
culation  of  the  production  rates  m,  and  hear  release  rate  H  will  be  introduced 
in  Sec.  4.1. 


2.2.2  Equations  of  State  for  Single  Species  and  Mixture 

The  compressible  Navier-Stokes  equations  (2.9-2.11)  are  open  equations 
since  the  number  of  unknown  variables  (p,  u,  E  and  P )  is  larger  than  the  num¬ 
ber  of  equations.  The  equation  of  state  (EOS)  is  needed  to  couple  the  pressure 
P  and  the  specific  energy  E  in  order  to  close  the  Navier-Stokes  equations.  In 
other  words,  we  require  functions 

P  —  P(E,p,u )  and  E  —  E(P,p,u).  (2-15) 

EOS  are  thermodynamic  equations  describing  the  state  of  matter  under 
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a  given  set  of  physical  conditions.  The  state  of  gas  is  commonly  described  by 
the  state  variables:  pressure  P ,  temperature  T  and  specific  volume  V  =  1/p. 
Actually,  experiments  show  that  only  two  state  variables  are  needed  to  define 
the  state  of  pure  gas  in  equilibrium  or  undergoing  a  steady  or  quasi-steady 
process.  In  other  words,  there  exists  a  function  f(P,  V,  T )  =  0. 

Specific  internal  energy  e  is  defined  as  the  energy  associated  with  the 
random,  disordered  motion  of  molecules.  The  e  is  a  function  of  the  state  or 
state  variables  (since  state  is  described  by  state  variables): 

e  =  e(V,T).  (2.16) 


Specific  enthalpy  H  is  defined  as 

H  =  e  +  PV  . 


(2.17) 


Using  the  first  law  of  thermal  dynamics,  it  can  be  derived  that  the  specific 
heat  for  a  constant  volume  process  is 


Cv  = 


(2.18) 


It  can  be  derived  that  the  specific  heat  for  a  constant  pressure  process  is 


a 


(2.19) 


The  Cp  and  Cv  are  thermodynamic  properties  and  depend  only  on  the  state. 
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Besides  the  internal  energy  e  and  the  enthalpy  H,  another  important 
thermodynamics  concept  is  the  entropy  S  defined  by  the  second  law  of  ther¬ 
modynamics.  Combining  the  first  and  second  law,  we  have  the  fundamental 
thermodynamics  identity: 

de  =  —PdV  +  TdS  +  T^UidNi  .  (2.20) 

Here  last  term  TluldNl  denotes  the  change  of  internal  due  to  changing  numbers 
of  particles.  Nt  is  the  number  of  particle  i  and  ut  is  the  chemical  potential  of 
partical  i. 


Our  model  scramjet  combustor  simulation  makes  an  assumption  of  ideal- 
gas  EOS.  The  ideal-gas  satisfies  the  law  of  Boyel  and  Gay-Lussac  (law  of  ideal 
gas): 


PV 


RT 

~w 


(2.21) 


where  R  is  the  universal  gas  constant  and  M  is  the  molecular  weight  of  the  gas. 
It  can  be  derived,  through  the  concept  of  entropy,  that  the  internal  energy  of 
an  ideal  gas  only  depends  on  the  temperature  of  the  gas  and  is  independent 
of  the  specific  volume,  i.e,  e  =  e(T).  The  enthalpy 


H  =  e(T)  +  PV 


e(T)  + 


RT 

~M 


(2.22) 


of  ideal-gas  is  also  a  function  of  the  temperature  only.  It  is  easy  to  see  Cp  and 
Cv  are  both  functions  of  T  only.  Since  internal  energy  e  only  depends  on  T, 
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we  have 


de 


[%)vdT+{w)TdV=c^dT 


since  the  second  term 


de 

av 


0.  We  can  also  have 


(2.23) 


dH  =  Cp(T)dT 


(2.24) 


Subtract  Eq.  (2.23)  from  Eq.  (2.24)  we  get 


dH  —  de  —  (CP(T)  -  Cv(T))dT 


(2.25) 


Since  dH  —  de  =  d(H  —  e)  —  d(PV )  =  d(^),  we  can  further  derive 


Cr(T)  -  CV(T)  =  A 


(2.26) 


The  ratio  of  specific  heats  is  denoted  by  7: 


7  = 


Cp 

Cv 


(2.27) 


The  specific  heat  CP(T )  of  all  ideal-gas  is  approximated  by  fourth  order  poly¬ 
nomials  with  the  coefficients  calculated  through  a  least  square  fit  by  NASA 
[43] 

__  _  rn  ,  rn 2  1  m3  i  rn4  23) 


— —  —  0\  +  02  T  +  O3T2  +  0.4T3  +  05  T4 
K 


H{T)  follows  the  integral  form  of  equation  (2.28) 


H 


a'2  m  ,  a3  m2  ,  rp3  .  a5  rp4  ,  °6 


- =  O!  +  —  T  +  —T  +  — +  —  T  +  — 

RT  2  3  4  5  T 
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(2.29) 


Here  Gq,  a 2,  03,  04,  05,  a6  and  07  are  the  numerical  coefficients  tabulated  in 
NASA  thermodynamic  table  [43].  The  pressure  P  and  the  internal  energy  E 
of  a  ideal  gas  can  be  coupled  by  equation: 

e  =  H(T)  +  PV  =  (ai  +  +  °^T2  +  °fT3  +  ^T4  +  ^  RT 

v  ’  \  2  3  4  5  T  M  J 

(2.30) 

1.  when  we  know  P  and  V,  we  first  calculated  the  temperature  as  T  — 
PVM/R ,  then  use  the  temperature  to  calculate  e. 

2.  when  we  know  e  and  V,  we  solve  the  non-linear  equation  of  e  =  H(T)  + 

for  the  temperature  of  T,  then  P  is  calculated  by  P  =  RT / (MV). 

The  specific  energy  E  is  related  to  the  specific  internal  energy  e  by  E  = 
e  +  p{u2  +  v2  +  w2)/ 2. 

The  mixture  of  ideal  gas  is  also  an  ideal  gas.  The  above  equations  can  be 
used  with  some  modifications.  To  adapt  the  EOS  of  pure  species  to  a  mixture, 
the  molecular  weight  M  is  replaced  by  the  average  M  for  mixture  of  n  species: 

n 

M  =  (2-31) 

i=  1 

The  specific  heat  Cp  and  enthalpy  PI  of  mixture  can  be  calculated  empirically 
by: 

n 

CP  =  J2  y‘cp,  •  <2-32) 

i—  1 

and 

n 

H  =  (2-33) 

i— 1 
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2.2.3  Transport  Properties 


The  transport  properties  included  in  the  our  model  scramjet  combustor 
are  viscosity,  thermal  conductivity,  and  mutual  diffusivity.  The  dynamic  vis¬ 
cosity  and  thermal  conductivity  are  calculated  by  semi-empirical  mixture  rules 
[10]: 


AO 

I  I  M,  \ "  YjGij  ’ 

^  Y  Mj 


_ h _ 

1  Mi  VJ=n  YiG*i  ’ 

^  Yi  Mj 


where  dimensionless  quantity  Gl3  are 


(2.34) 

(2.35) 


Gi 


1 

71 


-1/2 


(2.36) 


The  dynamic  viscosity  and  thermal  conductivity  of  a  single  gas  species  are 
computed  by  the  elementary  gas  model: 

Hi  =  2.6693  •  10-5^  .  (2.37) 

52P(2’2)  ' 

»  v 

ki  =  /I i  '  (2.38) 

Here  Mj  is  the  molecular  mass  of  species  i.  5,  is  a  characteristic  diameter  of 
the  molecule  i,  which  is  often  called  the  collision  diameter.  The  dimensionless 
quantity  is  called  the  “collision  integral  for  viscosity”  which  accounts  for 
the  paths  molecules  take  during  a  binary  collision.  Cp  is  the  heat  capacity,  the 
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calculation  of  which  has  been  introduced  in  Sec.  2.2.2. 


The  calculation  of  the  binary  diffusion  coefficient  of  mixture  is  complex 
and  is  omitted  here.  Details  of  the  binary  diffusion  coefficient  calculation  can 
be  found  in  [15]. 

2.3  Turbulent  Modeling  Approaches 

There  are  three  main  methodologies  for  turbulent  flow  simulation:  Reynolds- 
averaged  Navier-Stokes  (RANS),  large  eddy  simulation  (LES)  and  direct  nu¬ 
merical  simulation  (DNS).  All  these  three  methodologies  solve  the  governing 
Navier-Stokes  equations.  However  they  have  a  different  way  of  predicting  the 
fluid  dynamics. 

In  Sec.  2.1,  we  describe  three  different  length  scales  of  turbulence.  DNS 
is  the  simulation  methodology  which  resolves  the  turbulent  scale  down  to  the 
Kolmogorov  scales.  RANS  only  resolves  the  large  and  energetic  scales.  LES 
lies  midway  between  DNS  and  RANS.  The  grid  resolutions  of  these  three  kinds 
of  turbulence  simulation  approaches  are  illustrated  in  Figure  2.4.  We  introduce 
these  methodologies  and  discuss  their  advantages  and  disadvantages. 

2.3.1  DNS 

DNS  solves  the  Navier-Stokes  equations  directly  without  any  turbulence 
modeling.  In  DNS,  the  whole  range  of  spatial  and  temporal  scales  of  turbu¬ 
lence,  from  the  large  integral  scale  lc  down  to  the  small  Kolmogorov  scale  rjk, 
are  resolved  directly. 
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Direct  numerical  simulation  (DNS)  ADNS 

Resolved  Modeled 


Large  eddy  simulation  (LES)  t\LES 

Resolved^  Modeled  ^ 

Arans  Reynolds  averaged  Navier-Stokes  equations  (RANS) 

Figure  2.4:  Grid  resolution  of  DNS,  LES  and  RANS  [47] 

DNS  is  computationally  expensive.  We  know  from  Sec.  2.1  that  %  ~ 
.Re-3/4.  Since  DNS  needs  to  resolve  both  the  large  integral  scale  and  the  small 
Kolmogorov  scale,  the  number  of  grids  needed  to  resolve  the  small  scales  in 
each  dimension  is  ~  Re3/4.  Thus  the  total  number  of  grids  in  three  dimensional 
simulation  is  ~  Re9/4.  Assuming  the  time  resolution  is  proportional  to  space 
resolution,  the  overall  cost  is  about  Re3. 

This  huge  computational  cost  prevents  the  use  of  DNS  for  a  wide  variety 
of  flows,  especially  at  large  Reynolds  number.  For  a  three  dimensional  flow 
with  Re  ps  5000,  the  grid  points  needed  will  be  250  million  cells.  Thus  the 
application  of  DNS  is  restricted  to  flows  with  low  Reynolds  number. 

Although  DNS  is  expensive,  the  accuracy  of  the  DNS  result  can  be  com¬ 
parable  with  the  experiment.  In  circumstances  when  the  experiment  data  are 
unavailable,  the  DNS  data  is  often  used  for  validation  and  verification  purpose. 
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In  conclusion,  DNS  works  as  a  powerful  research  tool  for  turbulent  flows  at 
low  to  moderate  Reynolds  numbers. 


2.3.2  RANS 

The  RANS  solves  the  time-averaged  Navier-Stokes  equations.  In  RANS, 
the  fluid  dynamics  are  decomposed  into  two  parts:  the  time  averaged  part  and 
the  fluctuation  part.  There  are  two  kinds  of  decomposition:  the  Reynolds  aver¬ 
aging  decomposition  and  the  Favre  averaging  decomposition.  In  the  Reynolds 
averaging  decomposition,  an  instantaneous  solution  variable  0  is  decomposed 
into  time  averaging  part  0  and  fluctuation  part  (j) : 

<j)=(j)  +  (j)'  (2.39) 

In  the  Favre  averaging  decomposition,  an  instantaneous  solution  variable  0  is 
decomposed  into  mass  weighted  time  averaging  part  0  and  fluctuation  part 


0  =  0  +  0" 


(2.40) 


where 


(2.41) 


In  derivation  of  governing  RANS  equations,  Reynolds  averaging  is  applied 
to  density  p,  pressure  p ,  shear  stress  r  and  heat  flux  q\  while  Favre  averaging  is 
applied  to  the  velocity  held  u,  specific  internal  energy  e  ,  and  specific  enthalpy 
h.  The  Governing  Favre-Averaged  Navier-Stokes  Equation  for  compressible 
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flow  are: 


1.  continuity  equation 


dp  d 

at  +  te  (?“‘)  =  0 ' 


2.  momentum  equation 


m  (ik“}  +  £ +  (fy "  p“‘“0  ■ 


(2.42) 


(2.43) 


3.  energy  equation 


_0_ 

Ft 


p[  e  +  -UiUi  I  +  -pu'-u- 


+ 


d 

dxj 


P[h+  ^ UiUi  1  +  y  pw"«" 


Ui 


d_ 

dx 


3  L 


Ty  -  put  Uj  I  Ui  -  q3 


d 


dxi 


n  7  //  //  //  //  .  // 

-pUjh  -/»lillitli  ~P"iTL 


(2.44) 


Here  we  have  four  unknown  terms, 

1.  pu”u”,  Favre-averaged  reynolds  stress  tensor, 

2.  pu'-h" ,  turbulent  transport  of  heat, 

3.  pu”u”u ",  turbulent  transport  of  kinetic  energy, 

4.  pu'-Tij ,  turbulent  molecular  diffusion. 

These  terms  reflect  the  effect  of  turbulence  on  the  mean  flow. 

To  close  the  mean  flow  equations,  a  turbulence  model  is  needed  to  calcu¬ 
late  these  terms.  These  turbulence  models  are  often  based  on  the  Boussinesq 
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hypothesis  that  the  Reynolds  stress  could  be  linked  to  the  mean  rate  of  defor¬ 
mation: 

77  fau,  au,\ 

r =  +  (2.45) 

where  /.q  is  the  turbulent  viscosity. 

There  are  many  flavors  of  R.ANS  turbulence  model.  To  name  a  few,  there 
are  zero  equation  model,  the  mixing  length  model;  the  one  equation  models, 
like  the  Spalart-Almaras  model;  the  two  equation  models,  the  k  —  e  style 
models  and  the  k  —  u>  model  for  example;  the  seven  equation  model,  namely, 
the  Reynolds  stress  model.  The  number  of  equations  represent  the  number  of 
additional  PDEs  that  need  to  be  resolved.  More  details  of  these  R.ANS  model 
can  be  found  in  [47]. 

Compared  with  DNS,  R.ANS  can  obtain  an  averaged  solution  on  a  much 
coarser  grid.  Thus  R.ANS  is  very  economical  by  computational  resource  and 
time.  For  many  applications,  the  averaged  solution  or  the  steady  state  solution 
is  preferable.  Thus  R.ANS  is  a  widely  used  approach  for  engineering  problems. 

However,  R.ANS  could  not  predict  the  transient  behaviour  of  fluctuating 
flow  held.  In  problems  like  turbulent  mixing  and  combustion  where  the  chem¬ 
ical  reactions  are  driven  by  the  unsteady  turbulent  mixing  of  the  fuel  and  the 
oxidizer,  R.ANS  might  lose  its  accuracy. 

2.3.3  LES 

In  LES,  the  larger  turbulent  motions  are  resolved,  while  the  effects  of 
smaller  scale  motions  are  modeled.  LES  lies  in  the  midway  between  DNS 
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and  LES  in  terms  of  computational  expense  and  simulation  accuracy.  LES  is 
motivated  by  the  shortcomings  of  DNS  and  RANS.  Compared  with  RANS, 
LES  is  more  accurate  because  it  resolves  the  large  scale  unsteady  motions 
down  to  the  inertial  range.  Compared  with  DNS,  LES  is  more  economic  in 
computational  cost  because  DNS  spends  most  of  the  computational  resource 
in  resolving  the  small  dissipative  motions  while  LES  models  the  effect  of  these 
motions. 


The  governing  equations  of  LES  are  the  space  filtered  Navier-Stokes  equa¬ 
tion.  The  motions  and  eddies  with  scale  smaller  than  the  filter  width  are 
removed  by  the  filtering.  The  larger  and  important  eddies  remain  and  are 
governed  by  the  resulting  equations. 


A  space  filter  G  decomposes  a  variable  (j){x)  into  the  sum  of  a  filtered 
component  and  a  residual  component  <f> .  The  filtered  conponent  is  defined 
by 

*  =  (2.46) 

For  compressible  flow,  the  Favre  filtering  is  defined  in  the  similar  way  as  the 
Favre  Averaging  in  Sec.  2.3.2: 


J p(x  )4>{x  )G{x ,  x')dx 
f  p(x')G(x,  x')dx' 


(2.47) 


The  Favre  filtered  Navier-Stokes  equations  for  compressible  flow  are  [40]: 
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1.  Filtered  continuity  equation 


dp  dpUj 
d  t  dxj 


2.  Filtered  momentum  equation 


dpu.i 
dt  + 


dpUiUj 

dxj 


+ 


dP 

dxj 


d 

dxj 


3.  Filtered  energy  equation 


dt  dxj  dxj  dxj  \  dxj  I  dxj  dxj  I 


dqf  ,  ^ 


dxi 


+  Se 


4.  Filtered  concentration  equation 

dprfk  dpipkUj  _  _d_  f _  QqIS  ,  e 

dt  dxj  dxj  y  dxj  J  dxj  k 

The  subgrid  scaled  (SGS)  terms  TtJ,  q\H\  and  q ^  are  expressed  as 

Pj  =  p  (upPj  -  UjUj )  ,  (2.48) 

q\H)  =  P  (cpTui  -  CpTui'j  ,  (2.49) 

q^i  =  p  (?PkUi  -  •  (2.50) 

Here  the  filtered  variable  p,  ip,  ut  ,  p  and  E  denotes,  respectively,  the  fil¬ 
tered  density,  species  mass  fraction,  velocity,  pressure,  and  total  specific 
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onormr  P 


Hk 


P 

P 


(2.51) 


where  ek  is  the  hltered  specihc  internal  energy  of  each  species,  we  assume 
g-T')  =  0  [26].  The  viscous  stress  tensor  ftJ  is 


T, 


ij 


2duk  \ 
3dxkij) 


(2.52) 


where  fi  =  up  is  the  hltered  dynamic  viscosity. 

There  are  several  SGS  models  to  close  the  space  hltered  Navier-Stokes 
equation.  We  discuss  the  Classical  Smagorinsky  model  [62]  and  the  dy¬ 
namic  Smagorinsky  model  [24]  in  detail  in  this  section. 


Classical  Smagorinsky  Model 

The  Smagorinsky  model  assumes  a  the  local  equilibrium  of  the  subgrid 
scales.  The  Smagorinsky  model  for  eddy-viscosity  is  combined  with  the 
Yoshizawa  model  for  SGS  turbulent  kinetic  energy.  For  compressible 
how,  the  subgrid  terms  in  species  concentration  and  energy  equation 
also  use  the  gradient-transport  models. 

The  SGS  stress  is  decomposed  into  anisotropic  and  isotropic  tensors 
and  each  tensor  is  modeled  separately: 

6a  x  5u 

T»  (t  nr  v  \  i  t  LJ 

ij  \  ij  kk  g  )  i  -I-  kk  ^ 
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(2.53) 


with 


and 


an-isotropic  tensor 


isotropic  tensor 


rna  _  rri  _  rp  lj 

-L  in  in  -L  kk  ~ 


-lj 


rpi  _  rp  Uij 

J-nA  ~  J-kk.—T 


~*3 


(2.54) 


(2.55) 


The  anisotropic  part  is  modelled  based  on  the  concept  of  turbulent  vis¬ 
cosity 

ly  =  -2  «S-  (2.56) 

where  .S'.,  -  )  (yy  +  Tp)  ,  Sfj  =  Slf  -  Vt./  By  the  classical  Smagorin- 
sky  model,  the  simple  algebraic  model  holds  for  ut: 


ut  =  (C'.A)2^!  (2.57) 

with  A  =  ( AxAyAz )1/3  and  S  =  ^/2 SijSij.  Cs  is  a  constant  with  value 
between  0.065-0.25.  Thus  for  the  anisotropic  part  we  have 


Ty  =  -2p(C,A)2|S|S$  (2.58) 

The  isotropic  part  of  SGS  stress  based  on  Yoshizawa  model  [74]  with  the 
argument  that  ksos  =  (C7A)2|S|2, 

Tkk  =  -2p(CVA)2|£|2  (2.59) 
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Cs  is  a  model  constant.  The  full  SGS  terms  can  be  expressed  as: 


Ttj  =  -2p(C,A)2|S|S“  +  5-f2p(C,Kf\S\2  (2.60) 


A  gradient-transport  assumption  is  used  for  the  SGS  mass  species  trans¬ 
port, 


,  _  .  vt  dYk 

k y'  ^  Sct  dxj 


(2.61) 


where  Sct  is  the  turbulent  Schmidt  number.  The  SGS  heat  transport 
can  be  derived  with  similar  assumption, 


n.  =  _WidL  =  _jcdL 

'  Prt  dxj  1  dxj 


(2.62) 


where  Prt  is  the  turbulent  Prandtl  number. 


Once  the  four  coefficients  Cs,  Cj,  Prt  and  Sct  are  specified,  the  filtered 
Navier- Stokes  equation  for  compressible  flow  is  closed.  According  to  Ger- 
mano  et  al.  [24],  there  are  several  limitations  of  the  classical  Smagorin- 
sky  model.  First,  the  model  constant  Cs,  Ci ,  Prt  and  Sct  should  be 
changed  to  adapt  to  different  flows.  Second,  the  model  does  not  have  a 
correct  limiting  behaviour  neat  walls.  Third,  the  model  always  assumes 
a  positive  eddy  viscosity,  even  for  laminar  flows. 


To  overcome  these  limitations,  a  Dynamic  Smagorinsky  Model  (DSM)  is 
introduced  by  Germano  [24], 
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Dynamic  Smagorinsky  Model 


The  Dynamic  Smagorinsky  model  calculates  the  model  coefficients  Cs, 
Ci,  Prt  and  Sct  dynamically.  In  DSM,  a  dynamic  procedure  is  intro¬ 
duced  base  on  the  Smagorinsky  model  applied  at  two  different  filter 
levels.  The  SGS  stress  at  grid  level  A  is 


Tti  =  -2p(C,A)2|S|S“  +  ^2p(C/A)2|S|2  (2.63) 

A  test  filter  is  introduced  with  a  filter  width  A  larger  than  the  grid  filter. 
The  SGS  stress  at  the  test  filter  level  A  is  by  defination: 


tij  pUjUj 


pUipUi 


p 


(2.64) 


At  the  same  time  the  SGS  stress  at  the  test  filter  level  can  be  modeled 
by: 

Tti  =  -2?(C,A)2|f|f°  +  Sf2f(C, A)2|I|2  .  (2,65) 

By  applying  Germane* ’s  identity,  we  have 


L{j  t{j  Tij 


pUipUi 


pUipUi 


p 


p 


(2.66) 


where  Ltj  is  the  Leonard  stress  tensor.  The  anisotropic  part  of  is 


T  1  _  fa 

■Gj  Lij 


Tti  =  2(CsAYp\S\S^  -  2(CsA)2p\S\Sij  =  C'2SM«  (2.67) 
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where 


A2  (2.68) 

Eq.  (2.67)  includes  5  independent  equations  but  only  has  one  unknown 
variables  Cs.  Lily  [39]  proposed  a  least-squared  approach  to  solve  for  Cs : 


Ma-  = 

i] 


/5|5|5ft -2(A)  >'|S|S 


C: 


(MijM?,)  ’ 


(2.69) 


where  (.)  operator  denotes  the  average  over  equation  number.  If  the 
value  of  U2  is  negative,  its  value  will  be  clipped  to  zero. 


Actually,  the  value  of  C2  calculated  from  this  equation  is  quite  unstable 
with  space  and  time.  The  value  of  U2  may  be  zero  in  large  fraction  of  the 
flow  held  after  clipping.  To  remove  the  large  fluctuations  and  the  zero 
values  of  C2,  different  averaging  approaches  are  proposed.  Usually  the 
numerator  and  denominator  are  averaged  over  homogeneous  flow  direc¬ 
tions.  Meneveau  [46]  proposed  an  averaging  method  based  on  the  fluid 
imaginary  particle  trajectory.  Some  authors  used  a  time  averaging  of  C2 
which  might  introduce  errors  before  the  flow  becomes  statistical  station¬ 
ary.  In  the  absence  of  a  homogeneous  flow  direction,  local  averaging  may 
be  used.  In  our  model  scramjet  comustor  simulation,  the  local  averaging 
method  is  applied,  the  U2  is  clipped  to  be  zero. 


The  model  coefficient  Cj  is  computed  as 


C2  = 


(4J 

(M‘kk) 


(2.70) 
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where 


Mkk  = 


-p[2SijSij)+2[  |  )  p[2 SijSt3 


A2 


(2.71) 


Applying  similar  procedures  to  the  SGS  mass  species  transport  and  the 
SGS  heat  transport,  we  can  have  Prt  and  Sct  calculated  as: 


Prt 
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((Eif’AffV) 
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(2.72) 

(2.73) 

(2.74) 

(2.75) 

(2.76) 

(2.77) 


2.4  Numerical  Scheme 

In  our  model  scramjet  combustor  simulation,  a  conservative  finite  differ¬ 
ence  approach  is  used  as  the  numerical  discretization  scheme.  Since  the 


33 


governing  equations  in  our  project  are  the  conservative  of  mass,  momen¬ 
tum,  energy  and  weight  of  species,  the  conservative  methods  (include 
the  conservative  finite  difference  methods  and  the  finite  volume  meth¬ 
ods)  are  more  suitable  than  the  conventional  finite  difference  methods. 
The  conservative  methods  are  also  better  for  the  shock-capturing  prob¬ 
lems:  they  can  estimate  the  shock  speed  correctly  and  satisfy  the  the 
Rankine-Hugoniot  relation  across  discontinuities. 


In  conservative  finite  difference  methods,  the  governing  equations  can  be 
solved  dimension  by  dimension.  Thus  the  conservative  finite  difference 
methods  are  better  than  the  finite  volume  methods,  which  could  be  com¬ 
putationally  expensive  for  multi-dimensional  problems  with  high  order 
accuracy.  For  each  dimension,  we  have  the  equation 


dU  dF 
dt  dx 


(2.78) 


where  F  is  a  function  of  U .  The  discretization  is  simple: 


dUt  F,+j  -  F,-i 
dt  Ax 


0 


(2.79) 


where  Ut  is  the  cell-averaged  conservative  variable  of  cell  b  Fi±\  is  the 
flux  at  location  i  ±  |  (cell  faces  of  cell  i).  The  value  of  Fi±i  can  be 
approximate  F)’s  in  cell  centers  of  local  stencil  using  reconstruction  or 
interpolation. 


In  our  model  scramjet  combustor  simulation,  a  5th  order  finite  differ- 
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ence  Weighted  Essentially  Non-Oscillatory  (WENO)  scheme  of  Jiang 
and  Shu  [31]  is  used  as  our  conservative  finite  difference  scheme.  WENO 
schemes  use  the  idea  of  adaptive  stencils  to  automatically  adjust  the 
stencil  weight.  It  can  achieve  high  order  accuracy  and  non-oscillatory 
property  near  discontinuities  [31].  WENO  schemes  have  been  widely 
used  to  solve  compressible  Navier-Stokes  equations,  together  with  high- 
order  Runge-Kutta  methods. 

The  WENO  scheme  is  constructed  based  on  Essentially  Non-Oscillatory 
(ENO)  scheme.  Harten  [25]  constructed  the  first  ENO  scheme  in  1987. 
The  first  WENO  scheme  is  a  third  order  finite  volume  scheme  by  Liu, 
Osher  and  Chan  [41].  The  fifth  order  finite  difference  WENO  scheme 
applied  in  our  problem  is  constructed  by  Jiang  and  Shu  [31]  in  1996. 

The  reason  we  prefer  the  finite  difference  WENO  over  the  finite  volume 
WENO  is  that  the  finite  difference  WENO  is  more  economic  and  Fron- 
Tier  has  structured  grid.  Although  the  finite  volume  WENO  scheme 
can  allow  for  non-smooth  and  non-structured  scheme,  the  finite  volume 
WENO  scheme  are  several  times  expensive  than  the  finite  difference 
WENO  scheme  at  the  same  order  of  accuracy. 

The  ENO  scheme  applies  uniform  high  order  polynomial  reconstruction 

of  F- ,  i  on  local  stencils.  There  are  several  candidate  stencils  available 

,,=c2 

for  this  kind  of  polynomial  reconstruction.  ENO  chooses  the  smoothest 
stencil  and  use  that  stencil  to  reconstruct  flux  Fi+ i.  Near  discontinuity, 
ENO  choose  the  stencil  on  the  continuous  side.  Thus  the  ENO  scheme 
is  locally  adaptive. 
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Different  from  ENO,  WENO  makes  use  of  all  available  stencil.  Different 
weights  are  assigned  to  different  stencils  based  on  the  smoothness  of  each 
stencil.  In  smooth  region,  it  assigns  similar  weights  to  all  stencils  and 
could  achieve  a  higher  order  of  accuracy  than  ENO.  Near  discontinuity, 
WENO  assigns  small  weight  to  the  discontinues  stencils  and  large  weight 
to  continuos  stencil.  In  this  way,  the  WENO  scheme  is  also  locally 
adaptive. 

There  are  several  steps  in  of  WENO  reconstruction  [61].  Take  Eq.  (2.78) 
as  example. 

(a)  we  need  to  split  F(u)  into  two  monotone  parts: 

f(u)  =  f+  0)  +  f~  0)  (2.80) 

to  ensure  df+(u)/du  >  0  and  df~(u)/du  <  0.  The  Lax-Friedrichs 
splitting  is  commonly  used: 

f±(Ku)  =  ^(/(M)  ±au)  (2-81) 

where 

a  =  maxu\f  (w)|  (2.82) 

In  the  following  steps,  reconstruction  of  df+  ( u ) /du>  0  and  df~  (u) /du  < 
0  are  conducted  separately. 

(b)  The  approximation  of  Ft  from  each  sub-stencils  needs  to  be  com¬ 
puted.  Here  we  have  three  stencils  and  each  stencil  has  three  points. 
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The  approximations  are: 


?(0)  =-l  -If.  +  —  1 

i-\- 1/2  3^^  ^  1  '  g  i* 

(2.83) 

(1)  -  1  f  +  5  f  + 1  f 

i+1/2  6 ■'*-1  ^  6 3  Jt+1 

(2.84) 

fi+ 1/2  =  g/i  +  g/i+1  —  fi+2 

(2.85) 

(c)  The  weights  70,  71,  72  need  to  be  found  to  make  the  combination 

fi+ 1/2  =  To/i+1/2  +  T1//+1/2  +  72  4+1/2  (2-86) 

to  have  fifth  order  accuracy.  The  value  of  7’s  is  70  =  jq,  71  =  §, 
72  =  ^ 

(d)  Find  the  nonlinear  weights  u0,  07,  u2  to  make 


“*+1/2  ~  ^o/j+i/2  +  ^l/i+l/2  +  Ulfi+1/2  (2.87) 

fifth  order  in  smooth  regions  and  non-oscillatory  at  shocks.  The  a/s 
can  be  calculated  by 


UJk  = 


oik 


<+1  +  0?!  +  (+2 


(2.88) 


where 


Olk 


7  k 


(2.89) 


(e  +  ISk)P 

Here  e  is  introduced  here  to  prevent  zero  denominator.  Its  value  is 
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set  to  be  e  =  10  6. 


The  JS’s  are  calculated  by  taylor  expansion  analysis: 

13  1 

I  So  =  y^(/-2  —  2/_  i  +  /o)2  +  -(/-  2  —  4/_i  +  3/0)2 

(2.90) 

1 3  1 

ISi  =  ^(/- 1  -  2/0  +  A)2  +  i(/_1  -  h)2 

(2.91) 

13  1 

IS 2  =  Y^(/o  —  2/i  +  f2)2  +  -(3/o  —  4/i  +  f2y 

(2.92) 

The  time  discretization  of  WENO  schemes  is  realized  by  third-order 

TVD  Runge-Kntta  method: 

w(1)  =  un  +  A  tL(un) 

(2.93) 

m(2)  =  +  Iw(l)  +  IA  iL(w(!)) 

4  4  4 

(2.94) 

wn+1  =  +  ^A ,tL(uW) 

o  o  o 

(2.95) 

Thus  we  have  3rd-order  in  time  and  5th-order  in  space  finite  difference 
WENO  scheme. 

In  implementation  of  the  WENO  scheme  into  the  compressible  hyper¬ 
bolic  solver,  we  apply  eigenvalue  decomposition  to  the  hyperbolic  part  of 
Navier-Stokes  Equation  and  apply  WENO  scheme  to  decomposed  equa¬ 
tions  on  eigenvector  space. 

The  draw  back  of  WENO  is  that  it  can  lead  to  numerical  damping  of  tur¬ 
bulence,  since  these  fluctuations  can  be  seen  as  shock-like  oscillations  by 
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the  scheme.  In  many  scramjet  simulation,  authors  use  separate  schemes 
for  shock-capturing  and  the  non-shock  part.  ENO  or  WENO  is  used  near 
the  shock  while  other  schemes  are  used  at  other  place.  For  example,  in 
the  LES  simulation  of  hyshot-II[38],  Larsson  uses  unstructured  essen¬ 
tially  non-oscillatory  (ENO)  second-order  to  accurate  capture  shock  and 
a  HLLC  approximate  Riemann  solver  in  continous  region. 

2.5  Turbulent  Inflow  Generation 

When  LES  and  DNS  approach  are  used  for  turbulence  simulation  and 
the  turbulence  in  the  flow  field  is  mainly  resolved,  accurate  prescription 
of  inflow  turbulence  is  required.  There  are  many  effective  methods  to 
generate  inflow  turbulence.  These  methods  fall  into  two  categories:  syn¬ 
thetic  turbulence  techniques  and  rescaling-reintroducing  method.  The 
synthetic  turbulence  techniques  generate  a  random  field  using  models 
and  this  is  added  to  a  mean  flow  profile.  The  rescaling-reintroducing 
method  takes  the  flow  held  at  a  downstream  station,  rescales  it  and  im¬ 
pose  it  as  inflow  condition.  We  will  introduce  rescaling-reintroducing 
method  in  Sec.  3.5.2. 

The  random-noise  method  is  a  widely  used  synthetic  turbulence  tech¬ 
niques.  However,  the  turbulence  generated  by  the  random-noise  method 
lacks  energy  in  the  low  wavenumber  range.  As  the  result,  the  pseudo 
turbulence  would  damp  to  zero  quickly. 

Another  popular  synthetic  turbulence  technique  is  the  digital  filtering 
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approach  by  Klein  [33].  which  is  later  improved  by  Xie  [71]  and  Touber 
[65].  The  digital  filtering  method  of  Klein  [33]  assumes  that  the  second 
order  correlation  of  homogenuous  turbulence  has  a  Gaussian  form  while 
Xie  [71]  and  Touber  [65]  assume  it  has  a  exponential  form.  The  method 
use  here  follow  Touber’s  procedure. 


To  enforce  the  inflow  turbulence  to  have  a  second  order  correlation  of 
exponential  form,  the  digital  filtering  method  applies  the  discrete  filter 
operator  Fn 

N 

(2.96) 


Vk  =  FN(rk)  =  ^2  bjrk+j 

j=-N 


where  the  filter  coefficient 


h  ■ 

Uj  ~ 


bj 


2—/k=—N  uk 


(2.97) 


and 


bk  =  exp  (  -  — 'j 
n  J 


(2.98) 


The  new  set  of  variables  v  has  zero-mean  and  two-points  correlation  of: 


VkVk+m  (  nm 

=  exp 


VkVk 


2  n 


(2.99) 


Assume  vk  is  assigned  to  one  dimensional  flow  held  with  grid  size  of  Ax 
and  integral  length  scale  of  Ix.  We  set  n  to  be  Ix/Ax  and  N  to  be  2 n. 


40 


Then  the  correlation  function  R(xk  +  x)  of  one  dimensional  flow  field  is: 


R(xk  +  x)  —  R(xk  +  mAx )  =  h  k+l"  = 

VkVk 

(2.100) 

is  exponential  form.  Here  Xk  is  a  point  of  reference,  x  +  Xk  is  a  point 
some  distance  away  from  the  reference  point  and  m  =  x/ Ax.  Ix  is  the 
integral  length  scale. 

To  generate  a  two-dimensional  inflow  velocity  profile  for  three-dimensional 
flow,  the  two-dimensional  filter  coefficients  are  defined  as 

bjk  =  bjbk  •  (2.101) 

At  the  zero  time  step,  the  velocity  field  (fluctuation  part)  is  computed 
by  Eq.  (2.96).  At  the  next  step,  the  velocity  field  is  calculated  in  the 
same  way  and  then  correlated  in  time  to  the  previous  calculated  velocity 
field  by: 

vI;fc  =  v%ldexp  +  vkJl-exp  ■  (2.102) 

Here  is  the  velocity  field  for  the  new  step,  v^ld  is  the  velocity  field  of 
old  step  and  Vk  is  the  filtered  velocity  field,  At  is  the  time  step  and  r 
is  the  Lagrangian  time  scale.  Eq.  (2.102)  enforced  the  velocity  field  to 
follow  the  two-points  correlation  in  steam-wise  direction. 

For  turbulent  inflow  with  mean  inflow  velocity  of  u,  v,  w  and  Reynolds 


exp 


7r m 
2  n 


=  exp 


nx 

2L. 
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stress  value  of  R,  the  velocities  at  inflow  boundary  are: 


(2.103) 


where  uu,vv,ww  are  the  velocity  fluctuation  calculated  by  Eq.  (2.102). 

Experiment  shows  that  the  scramjet  model  combustor  has  an  inflow  tur¬ 
bulence  intensity  of  1.6%  (u  /U  where  u  is  the  rms  value  in  one  the  three 
components).  Considering  that  the  energy  in  the  randomized  hied  tends 
to  dissipate  very  quickly  and  we  are  using  WENO  solver,  a  turbulence 
inflow  with  density  of  2.3%  is  implemented.  We  have 

Ru  =  R22  =  R33  =  u'u  =  v  v  =  w  w'  =  1939m2/s2  (2.104) 

and 

Rij  =  0  (i  %  j)  (2.105) 

The  integral  length  scale  Ix  is  set  to  1.0  cm. 
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Chapter  3 


Turbulent  Boundary  Layer 


In  the  supersonic  ramjet  problem,  LES  is  used  to  simulate  the  turbulent 
flow.  Turbulence  gives  rise  to  eddies  of  different  length  scales.  The  large 
and  most  energetic  eddies  scale  with  the  bulk  velocity  and  the  integral 
length  scale  while  the  smallest  eddies  scale  with  the  dissipation  rate  and 
the  Kolmogorov  length  scale.  By  Kolmogorov’s  hypothesis,  these  small 
eddies  are  universal  and  can  be  modeled.  LES  computes  the  energetic 
and  flow  dependent  large  eddies  directly  and  models  the  small  eddies. 
LES  can  achieve  good  accuracy  with  a  lower  computational  cost  than 
Direct  Numerical  Simulation  (DNS)  at  high  Reynolds  numbers. 

LES  loses  its  power  for  flows  at  a  boundary  layer.  LES  assumes  a  sep¬ 
aration  between  the  large  energetic  eddies  and  the  small  eddies.  At  a 
boundary  layer,  the  large  and  energetic  eddies  scale  with  the  distance 
from  the  wall.  The  grid  size  that  resolves  the  main  stream  flow  is  not 
sufficient  to  resolve  such  energetic  eddies  at  the  boundary  layer.  To  re¬ 
solve  the  boundary  layer,  refined  grids  near  the  boundary  are  commonly 
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used.  Wall-resolved  LES  has  a  grid  scaled  with  the  small  eddies  at  the 
wall,  including  the  inner  most  part  of  the  boundary  layer.  However,  the 
computational  cost  to  resolve  the  inner  most  layer  can  be  as  expensive 
as  DNS.  Due  the  computational  burden  to  resolve  the  inner  layer,  it  is 
generally  accepted  that  the  inner  layer  should  be  modeled,  rather  than 
resolved. 

Despite  the  difficulty  of  resolving  the  boundary  layer,  the  ability  to  cap¬ 
ture  near  wall  processes  is  very  important  for  the  scramjet  simulation. 
Larsson  [36]  has  performed  a  Rayleigh-Fanno  analysis  of  the  Hyshot  com¬ 
bustor  at  nominal  conditions,  with  combustion  heat  release,  wall  friction 
and  wall  heat  losses  taken  from  RANS.  It  shows  that  the  pressure- change 
due  to  these  phenomena  is  about  +1  bar  (with  combustion  heat  release), 
+1  bar  (with  wall  friction),  and  -1  bar  (with  wall  heat  losses)  compared 
with  the  RANS  pressure  of  2.5  bar. 

3.1  Introduction  to  Turbulent  Boundary  Layer 

Turbulent  boundary  flow  is  different  from  a  laminar  boundary  flow.  In  a 
laminar  boundary  flow,  there  are  well-behaved  stream  lines.  In  a  turbu¬ 
lent  boundary  flow,  many  vortices  exist  in  the  near  wall  region.  Fig.  3.1 
shows  a  photo  of  a  fully  developed  Turbulent  Boundary  Layer  (TBL). 
The  existence  of  these  vortices  changes  the  mean  velocity  profile  at  the 
boundary.  Fig.  3.2  compares  the  velocity  profile  for  a  laminar  bound¬ 
ary  layer  and  a  turbulent  boundary  layer.  The  velocity  gradient  close 
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Figure  3.1:  Turbulent  boundary  layer  filled  with  eddies  of  many  different  scales 
(Image  source  J.  R„  Garcia  [22].) 

to  the  wall  is  much  steeper  for  a  TBL  than  for  a  laminar  boundary  layer. 


Figure  3.2:  Comparison  of  laminar  and  turbulent  velocity  profiles.  (Image 
source  Hoffman  [27].) 


In  this  section,  we  follow  the  exposition  of  Pope  [57].  A  TBL  has  several 
layers  and  each  layer  has  distinct  turbulent  characteristics.  These  layers 
are  shown  in  Fig.  3.3. 

Here  6  is  the  boundary  layer  thickness.  The  boundary  layer  thickness 
5  has  many  definitions.  The  most  widely  used  one  is  £99,  the  value  of 
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Figure  3.3:  Structure  of  a  turbulent  boundary  layer  [37]. 


wall-normal  distance  at  which  point  where  wall-parallel  velocity  reaches 
99%  of  the  free  stream  velocity. 


Close  to  the  wall,  the  kinetic  viscosity  v  and  the  wall  shear  stress  tw  are 
the  most  important  parameters.  Corresponding  viscous  velocity  scales 
and  length  scales  are  defined  based  on  the  the  kinetic  viscosity  v  and 
wall  shear  stress  rw.  The  friction  velocity: 

(3.1) 

is  the  viscous  velocity  scale.  The  velocity  measured  in  terms  of  the 
friction  velocity  is  denoted  by 


u 


u 

uT 


(3.2) 
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The  viscous  length  scale  near  the  wall  is  Sv  =  z//uT,  which  we  call  the 
wall  unit.  The  distance  from  the  wall  measured  in  wall-units  is  denoted 
by 


V  v 

t  =  —y  ■ 

Oi/  UT 


(3.3) 


y+  and  u+  are  dimensionless  variables.  y+  can  be  seen  as  a  local  Reynolds 
number.  It  determines  the  relative  important  effects  of  the  viscous  and 
turbulence  phenomena.  In  the  viscous  wall  region,  defined  as  y+  <  50, 
the  effects  of  molecular  viscosity  on  the  shear  stress  are  large;  while  in 
the  outer  layer,  defined  as  y+  >  50,  the  effects  of  molecular  viscosity  can 
be  neglected.  In  the  viscous  sublayer,  defined  as  y+  <  5,  the  turbulent 
shear  stress  is  much  smaller  than  molecular  shear  stress. 

According  to  Prandtl,  the  inner  layer  in  a  sufficiently  high  Reynolds 
number  flow  has  a  velocity  profile  determined  by  viscous  scale,  and  the 
inner  layer  profile  is  independent  of  the  main  stream  velocity.  In  other 
words,  u+  is  a  function  of  y+  only.  The  inner  layer  can  be  divided  into 
three  parts:  the  viscous  sublayer,  the  buffer  layer  and  the  log-law  layer. 
In  the  viscous  sublayer,  where  y+  <  5,  the  turbulent  shear  stress  can  be 
neglected  compared  to  the  molecular  shear  stress.  A  linear  relationship 
holds  for  the  viscous  sublayer: 


u+  =  y+  . 


(3.4) 


In  the  overlap  region  of  the  inner  and  outer  layers,  the  logarithmic  law 
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of  the  wall  of  Von  Karman  holds: 


u+  =  - ln(y+ )  +  B  .  (3.5) 

K 

Here,  B  is  a  constant  and  k  is  the  von  Karman  constant,  with  approxi¬ 
mate  values: 

k  —  0.41,  B  =  5.2.  (3.6) 

The  region  between  the  viscous  sublayer  and  the  log-law  region  is  called 
the  buffer  layer. 

3.2  Wall  Resolved  LES 

Wall  Resolved  LES  has  expensive  calculation  cost  for  flow  with  high 
Reynolds  number.  To  resolve  the  whole  boundary  layer  stucuture,  the 
Wall  Resolved  LES  need  to  resolve  the  dynamics  of  the  inner  layer  dom¬ 
inated  by  quasi-streamwise  vortices  which  have  sizes  in  same  order  with 
viscose  length  scale  5V.  To  resolve  the  inner  layer,  it  need  a  constant  grid 
spacing  scaled  with  wall  unit. 

Assume  that  the  domain  LES  resolves  in  the  inner  layer  has  a  size  of 
CiS  x  C2S  x  NSV  (We  assume  we  solve  the  inner  layer  from  the  wall  up 
to  the  position  of  y+  =  N).  To  resolve  the  inner  layer,  LES  should  have  a 
grid  size  of  A  ~  5V.  Thus  the  grid  needed  in  the  wall  parallel  direction  is 
C^/ A  ~  S/Sv  (i  —  1,  2)  cells.  According  to  Pope,  the  length  scale  ratio 
S/Sv  increases  approximately  as  0.09Re°'88.  Here  the  Reynolds  number 
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is  the  bulk  Reynolds  number  of  the  main  stream: 


v 


(3.7) 


In  the  wall  normal  direction,  we  need  N  cells.  Thus  the  total  number  of 
cells  needed  to  resolve  the  inner  layer  is 


AT 


(3.8) 


Thus  it  is  very  expensive  for  LES  to  resolve  the  inner  layer. 

In  contrast  to  the  inner  layer,  the  cost  to  resolve  the  outer  layer  is  much 
smaller.  In  the  inviscid  outer  layer,  the  energetic  motions  scale  with  the 
outer  length  scale  5.  Chapman  [13]  estimated  the  number  of  cells  need 
to  resolve  the  outer  layer  is  proportional  to  Re0A  for  flat  plat  boundary 
layer  and  independent  of  Re  for  plane  channel  flow. 

Fig.  3.4  shows  the  number  of  cells  needed  to  resolve  the  boundary  layer 
for  plane  channel  flow.  At  Reynolds  number  of  O(106),  about  99%  of  the 
points  are  used  to  resolve  the  inner  layer  (the  inner  10%  of  the  boundary 
layer).  For  high  Reynolds  numbers,  we  could  only  afford  the  coarse  grid 
size  which  only  resolve  the  out  layer  in  LES.  The  computational  cost  of 
the  coarse  grid  is  still  independent  or  weekly  dependent  of  Re. 

When  the  LES  is  implemented  on  coarse  grid  which  only  resolved  the 
outer  layer  and  does  not  resolve  the  inner  layer,  numerical  errors  rises  in 
the  near-wall  region.  As  we  seen  in  Fig.  3.2,  velocity  gradient  is  larger 
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Figure  3.4:  The  number  of  cells  needed  to  resolve  the  boundary  layer  for  plane 
channel  flow.  (Image  source  Piomelli  [53].) 
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near  the  wall.  Coarse  grid  LES  will  underestimate  the  wall  shear  stress 
and  distorts  the  exterior  LES  if  the  no  extra  model  is  applied  for  inner 
layer  LES.  Thus  wall  model  is  needed  for  the  inner  layer  on  the  coarse 
grid  simulation. 

3.3  Wall  Modeled  LES 

There  are  many  approaches  for  the  modeling  of  inner  layer  [53].  These 
models  are  divided  to  two  catergories:  wall  stress  models  that  calculate 
the  wall  shear  stress  and  hybrid  LES/RANS  method. 

Hybrid  LES/RANS  method  (left  plot  of  Fig.  3.5)  uses  the  RANS  to 
simulation  the  near-wall  region  and  LES  in  the  exterior  region.  The  grid 
resolution  tangential  to  the  wall  are  determined  by  the  exterior  LES,  the 
grid  resolution  normal  to  the  wall  are  determined  by  RANS.  The  hybrid 
LES/RANS  method  is  a  successful  and  widely  used  wall  model.  However, 
Hybrid  LES/RANS  method  is  numerically  more  expensive  than  wall- 
stress  model  because  full  evolution  equations  are  solved  down  to  the 
wall.  Another  drawback  of  hybrid  LES/RANS  method  is  that  the  skin 
friction  calculated  is  consistently  under-predicted  by  around  10%  —  15%. 

In  the  wall-stress  model,  the  LES  extends  all  the  way  down  to  the  wall. 
LES  only  resolves  the  outer  layer  motion.  The  inner  layer  motion  is 
resolved  by  a  simplified  wall  model.  The  wall  model  relate  the  instan¬ 
taneous  velocity  and  temperature  in  the  overlap  layer  from  the  LES  to 
the  wall  shear  stress  and  heat  flux,  and  provides  them  as  a  boundary 
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condition  back  to  the  LES. 


The  wall  model  can  be  algebraic  functions,  the  logarithmic  law  for  exam¬ 
ple,  or  differential  equations  (the  thin  boundary  layer  equation)  solved 
on  a  local  embedded  grid.  The  wall  model  used  in  our  model  scramjet 
combustor  simulation  is  the  equilibrium  wall  model  of  Larsson  [32], 


3.4  Equilibrium  Wall  Model 


The  equilibrium  wall  model  is  derived  from  the  thin  boundary  layer 
equation  with  the  assumption  that  the  change  of  physical  variables  along 
direction  tangential  to  the  wall  is  much  smaller  compared  with  the  change 
along  wall  normal  direction  and  we  can  introduce  the  approximation  that 
the  physical  variables  are  constant  along  wall  parallel  directions.  Another 
assumption  is  the  equilibrium  assumption  (the  pressure  is  constant).  As 
a  result  there  are  only  two  independent  variables  in  the  equilibrium  wall 
model  for  compressible  flow:  the  velocity  and  the  temperature.  The 
governing  equations  for  equilibrium  wall  model  are: 

=  0  ,  (3.9) 


d_ 

dy 


,  .  du 

+  '“)  Ty 


and 


d 


yt  y  \  d T  du 

P^,+Tr)l)y+{p  +  ,l‘)Udy 


=  0  . 


(3.10) 


dy 

Here  pw  and  rw  are  the  local  instantaneous  density  and  wall  stress.  The 
thin  boundary  layer  equation  and  the  equilibrium  wall  model  equations 
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are  simplified  RANS  equations.  See  Sec.  A  for  the  deviation. 

The  interaction  of  the  inner  layer  wall  model  and  outer  layer  LES  is 
shown  in  Fig.  3.5.  We  first  pick  a  point  in  the  overlap  region  of  the 
inner  layer  and  outer  layer  as  the  matching  point.  The  governing  one 
dimensional  equations  Eqs.  (3.9-3.10)  are  solved  on  embeded  grid  from 
the  matching  point  all  the  way  down  to  the  wall.  The  temperature, 
velocity  and  species  concentration  at  the  matching  point  from  LES  set 
the  outer  boundary  condition  of  Eqs.  (3.9-3.10).  The  no-slip  boundary 
condition  of  Eqs.  (3.9-3.10)  is  applied  at  the  Wall.  After  Eqs.  (3.9-3.10) 
are  solved,  we  obtain  the  tangential  velocity  and  temperature  gradient 
in  wall-normal  direction,  which  give  the  wall  shear  stress  and  heat  flux. 
Then  the  wall  shear  stress  and  heat  flux  are  fed  back  to  the  outer  layer 
LES. 

In  order  to  resolve  the  wall  normal  velocity  and  temperature  gradient, 
the  first  off-wall  point  should  have  a  small  distance  from  wall  in  the  same 
order  of  wall-units.  The  our  implementation,  the  first  node  is  at  y+  —  1. 
Then  geometric  grid  stretching  is  used  for  embedded  grid  from  wall  to  the 
matching  point.  In  this  way,  the  number  of  grid  points  N  from  the  wall 
to  the  matching  point  is  proportional  to  log  (Re).  It  is  weakly  dependent 
on  the  Reynolds  number.  In  our  simulation,  N  is  approximately  50. 

The  calculation  of  the  viscosity  coefficient  n  and  heat  capacity  cp  of  the 
gas  mixture  are  introduced  in  Sec.  2.2.3.  The  dynamic  viscosity  /q  is 
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Tw>i  TWji 


Figure  3.5:  Sketch  of  Wall  model  by  Larsson  [37].  Left  plot  shows  the  wal¬ 
l-stress  models.  The  filtered  Navier-Stokes  equation  are  solved  on  the  left  grid. 
The  wall  shear  stress  are  estimated  from  algebraic  relation  or  by  solving  thin 
boundary  layer  equation  on  embed  grid  (middle  grid).  Right  plot  demonstrate 
the  hybrid  LES/RANS  model.  The  Navier  Stokes  equations  are  solved  on  the 
right  grid  with  different  turbulence  models  for  LES  and  RANS  parts). 
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computed  by  the  mixing-length  model: 


/h 


—  exp 


(3.11) 


The  Von  Karman  constant  k  =  0.41,  and  A+  =  17.  Since  coefficients 
H,  Uti  and  cp  on  the  embedded  grids  are  dependent  of  the  solution  of 
one  dimensional  equations  Eqs.  (3.9-3.10),  an  iterative  solver  needs  to 
be  used.  We  start  from  an  initial  guess:  linear  profiles  of  u  and  T .  We 
compute  coefficients  /i,  /q  and  cp  on  each  point  based  on  the  initial  linear 
solution.  Then  we  solve  the  governing  equations  Eqs.  (3.9-3.10)  with  sin¬ 
gle  tridiagonal  matrix  algorithm  (TDMA)  sweep.  After  that,  we  obtain 
the  new  state  variable  u  and  T.  We  compute  the  new  coefficients  /i,  fit 
and  cp.  We  repeat  the  process  until  the  distribution  of  state  variables  u 
and  T  converge. 


After  we  have  converged  to  a  state  variable  distribution  u  and  T  in  the 
inner  layer,  we  compute  the  wall  shear  stress  and  heat  flux  by 


and 


7~w 


UO 

l^w  1 

yo 


(3.12) 


cp  T0  -  T,„ 

Qw  j-f  yw 

Pr  y0 


(3.13) 
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3.5  Grid  Sensitivity  Study  of  the  Turbulent  Bound¬ 
ary  Layer  Model 

3.5.1  Configuration 

In  this  section,  we  study  the  grid  sensitivity  of  the  boundary  layer  model 
for  2D  simulations  of  a  supersonic  flat  plate  boundary  layer.  The  inflow 
O2  has  an  average  density  of  1.241  xlO~4  kg/rn3,  velocity  of  1890  m/s 
and  temperature  of  1276  K,  same  with  the  scramjet  inflow  condition.  If 
we  have  laminar  inflow,  the  supersonic  flow  above  the  flat  plate  has  a 
transition  from  a  laminar  boundary  layer  to  a  turbulent  boundary  layer 
at  a  position  downstream  from  the  inlet.  The  transition  is  shown  in 
Fig.  3.6. 

However,  the  distance  from  the  inlet  to  the  transition  point  is  too  long  to 
be  computationally  affordable.  Thus  to  have  fully  developed  turbulence, 
the  rescaling-reintroducing  method  of  Urbin  and  Knight  [66]  is  used.  We 
discuss  this  method  in  Sec.  3.5.2. 

In  our  TBL  simulation,  the  wall  is  located  at  y  —  0.  The  computational 
domain  is  [240,  10]  mm  in  the  stream-wise  and  wall-normal  directions. 
The  recycling  station  is  located  at  40  mm  downstream  of  the  inflow.  To 
get  fully  developed  boundary  layer  turbulence  through  the  rescaling  and 
reintroducing  method,  the  simulation  was  run  for  about  20  flow  throughs. 
Here  5  is  the  inflow  boundary  layer  thickness  with  value  of  S  —  2.25  mm. 

The  error  in  the  calculation  of  the  wall  shear  stress  and  the  heat  flux 


56 


Laminar  flow 


Transition 


Fully  turbulent 


Figure  3.6:  Image  of  an  idealized  transition  process  on  a  flat  plate.  (Image 
source  Hoffman  [27].) 
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has  two  sources:  First,  the  model  error  of  the  turbulent  boundary  layer 
model;  Secondly,  the  discrepancy  between  the  real  state  variables  (u  and 
T)  and  the  numerical  value  of  these  values.  This  is  to  say  that  the 
physical  profile  of  the  matching  point  might  not  be  well  resolved. 

Larsson  [32]  argued  that  the  second  error  is  un-avoidable,  because  the 
first  grid  cell  off  the  wall  can  not  be  well  resolved.  Larsson  also  indicates 
that  we  could  greatly  reduced  the  second  error  by  not  using  the  first  grid 
cell  off  the  wall.  We  will  discuss  this  question  in  Section  3.5.3. 

3.5.2  Turbulent  Wall  Bounded  Flow  Generation 

We  have  discussed  the  digital  filtering  method  in  producing  the  free 
stream  turbulence  in  Sec.  2.5.  The  Rescaling-Reintroducing  method 
is  a  popular  method  to  introduce  inflow  boundary  turbulence.  The 
Rescaling- Reintroducing  method  was  first  proposed  by  Spalar  [63] .  Lund 
[42]  introduced  modifications  in  the  method  to  deal  with  the  growth 
terms.  Urbin  further  modified  the  method  to  adapt  better  to  the  com¬ 
pressible  turbulent  boundary  layer.  Here  we  follow  the  exposition  of 
Urbin  [66]. 

The  idea  of  the  Rescaling-Reintroducing  method  can  be  seen  in  Fig.  3.7. 
We  take  the  flow  field  at  a  downstream  station,  rescaled  it  and  impose 
it  as  inflow  condition.  In  the  rescaling  step,  flow  field  components  are 
decomposed  into  time  averaged  and  fluctuation  parts.  Separate  scaling 
laws  are  applied  to  each  part. 
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Figure  3.7:  Generation  of  inflow  boundary  condition.  The  plot  comes  from 
Urbin  [66] 

For  compressible  flow,  we  apply  the  Van  Driest-Fernholz  and  Finley 
transformation  of  the  velocity  U,  denoted  as  Uyd- 


(3.14) 


T,  U °°  ■  -if  A  U 

Uv  n  =  — —S'ln  A—— 


A 


Un 


and 


A  = 


(7  -  l)l2]MlPr, 


(3.15) 


(3.16) 


1  +  [(7  -  l)/2 ]MlPrt  ’ 
where  Prt  =  0.89.  The  inner  layer  and  outer  layer  are  scaled  in  different 
ways.  For  the  inner  layer,  we  have  by  the  law  of  the  wall: 


^vd  =  M^)/i(?/+)  • 


(3.17) 
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On  the  outer  layer,  we  have  by  the  defect  law: 


^VD  -  CvS  =  ,  (3.18) 

where  r/  =  y/5  is  the  outer  coordinate.  The  functions  f\  and  /2  are  both 
universal  and  are  independent  of  x.  Thus  we  have  the  scaling  law: 

tC.  i,„  =  puv nMyi) ,  (3.19) 

for  the  inner  region  and 

t/viD,  ini  =  PUv D,  rec (hinl)  +  (1  -  /Wd  ,  (3-20) 

in  the  outer  layer.  Here  the  subscript  ini  denotes  the  inlet  and  rec  denotes 
the  recycling  station.  /3  is  the  ratio  of  friction  velocity  at  the  inlet  station 
and  the  velocity  at  the  recycled  station.  The  treatment  of  the  fluctuation 
part  is  similar.  The  velocity  fluctuation  in  the  inner  and  outer  sides  are 
rescaled  by: 

=  Pu'L  (2/m>  *)>  U"T  =  /^recChinl,  t)  .  (3.21) 

The  rescaling  of  the  mean  wall-normal  velocity  and  the  temperature  is: 

K7  =  VSF  =  ^ec(hini)  ,  (3.22) 


60 


and 


T^  =  Tve  c(y+),  T^  =  Tie M. 


(3.23) 


Here  we  assume  that  the  steamwise  pressure  gradient  is  small  compared 
with  the  wall-normal  temperature  gradient.  The  fluctuation  td™1  of  the 
wall-normal  velocity  is  calculated  in  the  similar  way  of  u".  Assuming 
the  pressure  fluctuations  are  assumed  to  be  small  compared  with  tem¬ 
perature  fluctuations,  the  fluctuation  of  temperature  is  scaled  as: 


Thr  =  Trec(y+vz,t), 


Tinl^  ~  ^rec(77inl)  Zi  t) 


(3.24) 


Finally,  we  take  a  weighted  average  of  the  inner  and  outer  profile  for  the 
complete  profile  of  velocity  and  temperature: 


=  (C7  +  «inD  [1  -  W(r, ini)]  +  (E/£ T  +  <T)W(n inl)  ,  (3.25) 


with  the  weight 


W(rj) 


1 

2 


4(77  -  B) 

(1  -  2  B)rj  +  B 


/tanh{  4) 


(3.26) 


3.5.3  Choice  of  the  Matching  Point 

In  the  calculation  of  the  wall  shear  stress  and  heat  flux,  the  first  grid 
cell  off  the  wall  is  commonly  used  as  the  matching  point.  Larsson  [32] 
demonstrated  that  there  are  persistent  numerical  and  subgrid  errors  in 
the  matching  point  if  we  use  this  choice,  while  a  better  choice  can  elim- 
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inate  this  persistent  error. 


Assume  the  wall  is  located  at  y  =  0  and  the  matching  point  is  located  at 
y  =  ym-  Assume  ym  is  located  at  y+  >  50,  in  the  log  layer.  The  length 
scale  of  motions  in  the  log  layer  is  proportional  to  the  wall  distance  y. 
Assuming  the  stress  carrying  motion  has  length  scale  Lj  =  Cty  (i  =  0, 
1,  2)  and  assuming  that  we  need  N  cells  located  within  this  length  scale 
to  resolve  the  stress  carrying  motion,  we  have 


A Xi  < 


Li 

N 


(3.27) 


If  we  use  the  first  grid  off  the  wall  as  the  matching  point,  which  means 
ym  =  Ay,  we  will  have 

(3.28) 

which  means  C\  >  N.  However,  the  numerical  Nyquist  criterion  indi¬ 
cates  N  >  2.  And  the  kinematic  damping  by  the  wall  indicates  C\  <  2. 
Thus  Eq.  (3.28)  does  not  hold.  It  means  that  the  dynamics  at  the  first 
grid  off  the  wall  could  not  be  well  resolved.  The  TBL  model  is  fed  with 
inaccurate  information  if  we  choose  the  first  grid  cell  off  the  wall  as  the 
matching  point. 

Larsson  observed  that  it  is  not  required  that  the  TBL  to  be  applied 
between  the  first  grid  point  and  the  wall;  actually  the  TBL  model  equa¬ 
tions  are  valid  for  any  interval  from  the  wall  to  a  point  in  the  inner  layer. 
Larsson  found  that  he  could  achieve  grid  convergence  by  fixing  ym  and 
refining  the  grid. 
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3.5.4  Results  and  Analysis 


To  study  the  grid  convergence  of  the  wall  shear  stress  and  the  heat  flux 
as  the  resolution  is  varied,  we  conduct  numerical  experiments  with  four 
levels  of  resolution,  see  Table  3.1.  At  each  level  of  grid  resolution,  we 

Table  3.1:  Grid  resolution  study  of  the  turbulent  boundary  layer 


simulation 

Ax  (mm) 

Ay  (mm) 

grid  level  I 

0.5 

0.5 

grid  level  II 

0.25 

0.25 

grid  level  III 

0.125 

0.125 

grid  level  IV 

0.125 

0.0625 

conduct  two  simulations  with  two  choices  of  the  matching  point  location: 
the  first  choice  is  to  set  the  matching  at  a  fixed  location  of  ym  =  0.25 
mm  while  in  the  second  choice  the  matching  point  is  the  first  cell  off  the 
wall. 

The  average  shear  stress  and  heat  flux  at  the  wall  in  each  simulation 
are  shown  in  Fig.  3.8.  The  shear  stress  and  heat  flux  are  averaged  from 
the  recycling  station  to  the  outlet.  The  averaged  boundary  thickness  is 
about  3.0  mm  from  the  recycling  station  to  the  outlet,  corresponding 
to  a  averaged  inner  layer  thickness  of  about  0.3  mm.  Thus  the  fixed 
matching  poing  ym  =  0.25  mm  is  located  in  the  inner  layer. 

The  averaged  wall  shear  stress  and  heat  flux  converge  with  grid  refine¬ 
ment  when  the  location  of  the  matching  point  is  fixed  (y  =  0.25  mm). 
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(a)  Averaged  shear  stress  at  the  wall.  The  dashed  dot  horizontal 
lines  indicate  a  ±10  %  variation  about  the  averaged  shear  stress 


of  the  finest  grid. 


(b)  Averaged  heat  flux  at  the  wall.  The  dashed  dot  horizontal 
lines  indicate  a  ±10  %  variatio§4ibout  the  averaged  heat  flux  of 
the  finest  grid. 


The  averaged  wall  shear  stress  and  heat  flux  computed  from  use  of  the 
first  grid  cell  off  the  wall  as  matching  point  do  not  converge.  This  result 
verifies  the  point  made  in  Sec. 3. 5. 3:  the  velocity  and  temperature  pro¬ 
files  at  first  grid  cell  off  the  wall  are  not  well  resolved  and  thus  are  not 
accurate;  by  having  the  matching  point  several  cells  away  from  the  wall, 
we  could  achieve  convergence  of  wall  shear  stress  and  heat  flux. 

The  average  velocity  profiles  (Fig.  3.9)  of  four  grid  resolutions  with  the 
matching  point  fixed  at  y  =  0.25  mm  also  show  the  convergence  trend. 
The  velocity  profiles  of  grid  1 1- IV  agree  quite  well  with  each  other. 


Table  3.2:  Error  ratio  of  the  wall  shear  stress  and  heat  flux  of  each  simulation. 
ym  denotes  the  location  of  the  matching  point.  By  the  simulation  result  of  this 
section,  we  have  an  error  approximation  of  the  TBL  model  at  different  grids. 


Ay  (mm) 

Vm  (mm) 

number  of  cells  from 

number  of  cells  in 

error  ratio  in 

error  ratio  in 

ym  to  the  wall 

boundary  layer 

wall  shear  stress 

wall  heat  flux 

0.5000 

0.2500 

1 

6 

0.189 

0.133 

0.2500 

0.2500 

2 

12 

0.104 

0.028 

0.1250 

0.2500 

3 

24 

0.004 

0.046 

0.0625 

0.2500 

5 

48 

0.5000 

0.2500 

1 

6 

0.189 

0.133 

0.2500 

0.1250 

1 

12 

0.142 

0.053 

0.1250 

0.0625 

1 

24 

0.111 

0.008 

0.0625 

0.0312 

1 

48 

0.194 

0.148 

The  error  ratio  of  the  wall  shear  stress  and  heat  flux  of  each  simulation 
are  summerized  in  Table  3.2.  Here  the  marginal  convergence  of  the  av¬ 
eraged  wall  shear  stress  and  heat  flux  are  accessed.  Clear  convergence 
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(a)  Fixed  matching  point  at  y  =  0.25  mm 


(b)  Use  first  cell  off  the  wall  as  matching  point 


Figure  3.9:  Average  velocity  profiles  for  different  grids  in  the  boundary  layer. 
Uinf  is  the  free  stream  velocity.  Note  sfffecifically  the  compare  of  Grid  II  to  III 
in  the  range  of  0.2  <  y/S  <  0.6.  The  importance  of  (a)  over  (b)  is  subtle  but 


might  require  one  more  level  of  mesh  refinement.  We  observe  from  Ta¬ 
ble  3.2  : 

(a)  At  grid  level  IV  and  grid  level  III,  when  grid  size  Ay  is  less  or  equal 
to  0.125  mm  and  we  have  >  3  cells  from  the  matching  point  to  the 
wall  the  ratio  of  discrepancy  is  about  5%  for  both  the  wall  shear 
stress  and  the  heat  flux.  The  average  wall  shear  stress  and  heat 
flux  have  achieved  marginal  convergence. 

(b)  At  grid  level  II  when  grid  size  Ay  is  less  or  equal  to  0.25  mm  and 
we  have  2  cells  from  the  matching  point  to  the  wall,  the  ratio  of 
error  is  about  10%  for  both  the  wall  shear  stress  and  the  heat  flux; 

(c)  Using  the  first  cell  of  the  wall  as  the  matching  point  introduces 
an  consistent  error  of  about  20%  in  the  averaged  wall  shear  stress 
occurs  even  at  the  finest  grid  level. 

As  a  conclusion,  we  obtain  a  decreased  error  in  the  calculation  of  wall 
shear  stress  and  heat  flux  as  we  have  more  grid  resolution  within  the 
boundary  layer  and  more  cells  located  between  the  wall  and  the  match¬ 
ing  point  wall  normal  direction.  To  obtain  a  better  accuracy  in  the 
calculation  of  the  wall  shear  stress  and  heat  flux,  it  is  better  to  have 
local  mesh  refinement  at  the  boundary.  If  a  maximum  of  10%  error  is 
allowed  in  the  averaged  values  of  wall  shear  stress  and  heat  flux,  it  is 
predicted  that  a  mesh  resolution  of  Ay  >  0.25  mm  is  required  for  the 
scramjet  simulation. 

On  the  choice  of  the  matching  point,  we  should  abandon  the  use  of  the 
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first  grid  off  the  wall.  The  only  requirement  for  he  matching  point  is  that 
it  should  be  located  in  the  inner  layer,  the  matching  point  could  located 
from  the  first  cell  up  to  the  nth  cell  off  the  wall  where  n  =  0.15/ A. 
Larsson  estimated  using  the  fourth  cell  off  the  wall  is  sufficient  to  have 
accurate  input  at  the  matching  point.  Thus  we  can  choose  the  matching 
point  location  to  be 


n  =  mm 


(3.29) 


Chapter  4 


Turbulent  Combustion 


In  gaseous  combustion,  two  classes  of  flames  are  often  considered:  pre¬ 
mixed  and  non-premixed  (diffusion)  flames.  A  premixed  flame  is  a  flame 
in  which  the  oxidizer  has  been  mixed  with  the  fuel  before  it  reaches  the 
flame  front.  For  example,  combustion  in  homogeneous  charge  spark  ig¬ 
nition  engines  are  premixed  flames.  A  diffusion  flame  is  a  flame  in  which 
the  oxidizer  combines  with  the  fuel  by  diffusion.  In  diffusion  flame, 
reactants  are  initially  separated,  and  reaction  occurs  only  at  interface 
between  fuel  and  oxidizer.  One  example  of  a  diffusion  flame  is  a  candle 
flame.  The  flame  in  scramjet  combustor  is  also  a  diffusion  flame. 

In  this  chapter,  we  study  the  grid  sensitivity  of  our  finite  rate  chemistry 
model  in  the  context  of  one-dimensional  and  two-dimensional  flows.  The 
grid  sensitivity  analysis  in  three  dimensional  simulation  is  too  expensive 
for  a  detailed  exploration  of  all  its  many  parameters.  Based  on  ID  and 
2D  simulations,  we  excluded  some  model  and  parameter  choices.  In  ID, 
we  study  the  grid  sensitivity  in  resolving  the  premixed  and  non-premixed 
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laminar  flames  with  detailed  and  reduced  mechanism  separately.  In  2D, 
we  study  the  the  gird  sensitivity  of  turbulent  combustion.  Then  we  use 
the  turbulent  diffusion  coefficient  computed  in  a  2D  and  3D  context  in 
the  ID  problem  analysis  to  and  get  a  less  stringent  mesh  requirement. 


4.1  Chemical  Kinetics 


In  balance  equations  for  the  mass  fraction  of  species  i: 


d_ 

dt 


+ 


d 

dxj 


dx2 


+  m  , 


(4.1) 


2 1  yP^ij  is  the  local  rate  of  change  and  yfnijYt j  is  the  convection 

term.  ^pdjY^J  is  the  diffusive  flux.  rht  is  the  chemical  source  term, 
the  calculation  of  which  involves  the  chemical  Kinetics. 


The  reaction  rate  Wj  for  single  reaction  j  follows  the  Law  of  Mass  Action: 
the  reaction  rate  is  proportional  to  the  products  of  the  concentration 
reactants.  The  reaction  rate  Wj  for  reaction  j  in  a  mechanism  with  n 
reactions  is 


// 


Here  kfj  and  krj  represent  forward  and  reverse  rate  coefficients  of  reaction 
j  as  a  function  of  the  temperature,  and  are  the  stoichiometric 
coefficients  of  species  i  for  the  reaction  and  production  side  of  reaction 
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j.  The  forward  and  reverse  rate  coefficients  kfj  and  krj  can  be  computed 
empirically  by  the  Arrhenius  Law.  The  forward  rate  coefficient  is: 


%  =  A,r‘.  exp 

where  Aj,  rij  and  Ej  are  the  pre-exponential  constant,  the  temperature 
exponent,  and  the  activation  energy  of  reaction  j.  Their  values  can  be 
found  in  the  reaction  mechanism  tables.  The  reverse  rate  coefficient  krj 
is  computed  from  the  forward  coefficient  and  the  equilibrium  constant 
from  reaction  j  by 


(4.3) 


kfj 


n 


(4.4) 


where  Pa  =  1  bar.  The  A Hj  and  A Sj  are  respectively  the  enthalpy  and 
entropy  changes  of  reaction  j. 

The  chemical  source  term,  the  production  rates  rhi,  is  the  sum  over  all 
production  terms  in  the  mechanism, 

N 

m  =  ^  MiVjiWj  ,  (4.5) 

j=  i 

where  N  is  the  number  of  reactions,  AL,  is  the  molecular  mass  of  species 
i,  Wj  is  the  rate  for  reaction  j  and  zq-j  is  the  stoichiometric  coefficient  of 
species  i  in  the  reaction  j. 

The  heat  release  rate  is  expressed  as  the  sum  over  the  heat  release  of  all 
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chemical  reactions: 


N 

Hi  =  MihiWj  (4.6) 

j= 1 

where  ht  is  the  heat  release  in  reaction  i. 

In  last  several  decades,  lots  of  experiment  has  been  done  to  measure 
the  rate  of  gas  phase  reactions.  Many  authors  have  selected  reactions 
and  their  rate  data  and  combined  them  elementary  reaction  tables.  In 
our  problem,  the  fuel  is  H2.  There  are  many  studies  of  the  detailed 
hydrogen-air  reaction  mechanisms,  such  as  Jachimowski’s  33-reaction 
mechanism  [30],  Hong  et  al.’s  21-reaction  mechanism  [28]  and  Williams’s 
21-reaction  mechanism  [70].  Additional  discussion  can  be  found  from 
the  text  and  references  cited  in  [2,  5,  73,  59,  16].  The  mechanism  by 
[28]  has  8  species  (H,  O,  OH,  H20,  H2,  02,  H02,  and  H202)  and  20 
reactions,  as  in  Tab.  4.1.  The  rates  of  each  reaction  with  H2  and  02 
mixed  stoichiometrically  under  temperatures  of  1276  K  and  a  pressure 
of  100  kPa  are  shown  in  Fig.  4.1. 

In  the  detailed  mechanism,  the  chemical  source  terms  Wi  s  contain  the 
contributions  from  many  fast  reactions.  Thus,  the  reactive  chemistry 
equations  contain  a  system  of  stiff  non-linear  equations.  Solving  these 
non-linear  equation  requires  an  adaptive  solver  and  high  computational 
cost  [60]. 

Therefore,  the  reduced  chemistry  mechanisms  are  developed  and  are 
widely  used.  Balakrishnan  [3]  conducted  a  numerical  investigation  of 
the  extinction  and  ignition  limits  in  laminar  non-premixed  counter  flow 
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(a)  full  range  of  all  reactions 


(b)  enlargement  showing  minor  reactions 


Figure  4.1:  Reaction  rates  for  each  reaction  in  the  full  chemistry  mechanism 
for  H2/O2  mixed  stoichiometrically  under  model  scramjet  combustor  temper¬ 
atures  of  1276  K  and  pressure  of  100  kPa. 


for  both  detailed  and  reduced  chemistry  of  H2/Air  reactions.  Pantano 
[48]  used  a  four-step  reduced  mechanism  in  a  3D  direct  numerical  simu¬ 
lation  (DNS)  of  a  spatially  evolving  planar  turbulent  reacting  jet  for  the 
combustion  of  methane  with  air. 

The  reduced  reaction  mechanism  assumes  the  Quasi-Steady-State  for 
intermediate  species.  The  fast  reactions  depleting  the  quasi-steady  state 
intermediate  species  are  eliminated.  The  slow  rate  reactions  remain  and 
determine  the  rates  of  global  reactions. 

A  four-step  reaction  mechanism  [8]  for  H2/O2  interactions  has  been  suc¬ 
cessfully  verified  against  the  result  of  detailed  chemistry  in  the  main 
reaction  zone  for  high-temperature  ignition  [60] .  This  four-step  reaction 
mechanism  has  15  elementary  steps  and  the  4  global  steps. 
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I)  H  +  02  ^  OH  +  O 

II)  O  +  H2  +  M  ^  H  +  OH  +  M 

III)  OH  +  H2  ^  H  +  H20 

IV)  H  +  H  +  M  ^  H2  +  M 
with  global  reaction  rates  as: 

W/  =  Wi  +  we  +  W12  +  W15  -  W17  +  W 18  -  W20  +  W21 

wn  =w2-W3  +  we  +  W 12  +  W 14  +  wle  +  Wn  +  W 18  +  W21 

Win  =  W3  +  W4,  +  W8  +  W10  +  W13  +  W15  +  Wig 

Wjv  =  w5  +  w9  +  W10  -  W12  +  w16  +  W 17  +  W20 

It  neglects  H202  and  compute  the  partial  density  of  H02  from  quasi¬ 
steady  assumption.  It  aggregates  the  reaction  rates  of  elementary  re¬ 
actions  into  4  global  reactions.  This  four-step  reaction  mechanism  is 
used  in  the  finite  rate  chemistry  model  in  our  model  scramjet  combustor 
simulation. 

4.2  Flame  Structures 

4.2.1  Premixed  Flame 

For  a  premixed  flame,  the  fuel  and  oxidizer  are  completely  mixed  be¬ 
fore  the  combustion  takes  place.  The  fuel  and  oxidizer  are  mixed  at  low 
temperature  before  they  enter  the  combustion  chamber.  The  chemical 
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reaction  chain  is  sensitive  to  the  temperature.  At  low  temperature,  the 
“chain-breaking”  mechanism  drives  the  reaction  chain  thus  the  combus¬ 
tion  reactions  are  “frozin”.  When  a  strong  heat  source  (ignitor)  raises 
the  temperature  beyond  the  ignition  temperature,  the  combustion  starts. 

The  equivalence  ratio  of  a  given  mixture  is  an  important  parameter  for 
premixed  gases.  It  is  defined  as 


=  s 


Yf 

Y( 


(4.7) 


o 


where 


s  = 


v0M0 


(4.8) 


vFMp 

is  the  mass  stoichiometric  ratio.  Here  the  index  F  and  O  stands  for  fuel 
and  oxidizer  respectively.  Rich  combustion  happens  when  0  >  1  (the 
fuel  is  in  excess)  and  lean  combustion  happens  when  0  <  1  (the  oxidizer 
is  in  excess). 


In  premixed  combustion,  there  are  two  stable  states,  the  unburnt  states 
and  the  burned  gas  states.  They  are  separated  by  the  flame  front  as 
in  Fig  4.2.  The  laminar  burning  velocity  sF)  defined  as  the  velocity  at 
which  flame  front  propagate  in  the  direction  normal  to  itself  and  relative 
to  the  flow  into  the  unburnt  mixture  [51]. 

A  basic  case  of  a  premixed  flame  is  a  one-dimensional  laminar  premixed 
flame.  Computing  one-dimensional  laminar  premixed  flames  can  be  seen 
as  the  first  step  toward  more  complex  flames.  The  structure  of  one 
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unburnt 

mixture 


Figure  4.2:  Premixed  flame  (Image  source  [51]). 
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T,  Yj,  or  (bT 


(a)  Structure  of  a  laminar  premixed  flame(b)  Structure  of  a  laminar  diffusion  flame.  [1] 

[35]. 


Figure  4.3:  General  Flame  structure. 


dimensional  laminar  premixed  flame  is  shown  in  Fig.  4.3(a). 


4.2.2  Diffusion  Flame 

In  a  diffusion  combustion,  the  fuel  and  oxidizer  enter  combustion  cham¬ 
ber  separately  before  they  mix  and  burn  during  continuous  inter-diffusion. 
The  combustion  occurs  at  the  interface  between  the  fuel  and  the  oxidizer. 
The  structure  of  a  diffusion  flame  is  shown  in  Fig  4.3(b).  In  diffusion 
flames,  the  simple,  measurable  parameters,  like  the  rate  of  burning  and 
the  flame  velocity  in  premixed  flames,  can  not  be  easily  defined.  The 
burning  in  a  diffusion  flame  depends  more  on  rate  of  mixing  than  on 
the  rates  of  chemical  processes  involved.  The  rate  of  reaction  is  directly 
related  to  the  amounts  of  fuel  and  oxidant  diffusing  into  the  reaction 
zone. 
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4.3  Combustion  Models 


There  are  several  turbulent  combustion  models  for  the  premixed  and 
the  diffusion  flames.  For  the  premixed  flames,  there  are  Bray-Moss- 
Libby  Model  and  Coherent  Flame  Model  which  assume  the  infinitely  fast 
chemistry;  there  is  the  flamelet  model  based  on  G-eqnation  with  finite 
rate  chemistry  assumption.  For  diffusion  flames,  there  are  conserved 
scalar  equilibrium  model  (with  infinitely  fast  chemistry  assumption),  the 
flamelet  model  base  on  mixture  fraction  and  conditional  moment  closure 
model  (with  finite  rate  chemistry  assumption).  The  PDF  Transport 
Equation  model  and  the  linear  Eddy  Model  work  for  both  diffusion  flame 
and  premixed  flame.  We  will  discuss  the  flamelet  model  for  diffusion 
flames  in  detail.  The  description  of  other  models  can  be  found  in  [51]. 

William  [68]  saw  the  turbulent  diffusion  flames  as  an  combination  of 
stretched  laminar  flamelet.  Peter  [50]  derived  the  flamelet  equation 
Eq.  (4.10)  for  the  diffusion  flame.  The  flamelet  is  a  thin  reactive-diffusion 
layer  in  the  turbulence  flow  field  [51].  The  flame  surface  is  defined  as  the 
surface  of  stoichiometric  mixture  where  [51]: 

Z(x,t)  =  Zst  (4.9) 

Here  Z  is  the  mixture  fraction.  The  flamelet  equation  describe  the 
reactive-diffusive  structure  at  the  vicinity  of  the  flame  surface.  It  as¬ 
sumes  the  mass  fraction  and  temperature  can  be  expressed  as  a  function 
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of  the  mixture  fraction  Z.  The  governing  flamelet  equation  is: 


dip  p  x  d2ip 

P  dt  Le  2  dZ 2 


(4.10) 


where  ip  is  the  “reactive  scaler”  [51]  which  could  be  the  mass  fraction 
of  each  chemical  species  or  the  temperature.  If  ip  stands  for  the  mass 
fraction  of  species  i,  D  is  the  mass  diffusivity  of  the  species  i,  Le  is  the 
Lewis  number  Le,  =  A /{pcpDi)  =  D/Di ,  and  uj  is  the  chemical  produc¬ 
tion  rate.  If  0  represent  the  temperature,  D  is  the  thermal  diffusivity, 
Le  is  1  and  uj  is  the  heat  release  rate  in  chemical  reactions,  x  is  the 
scaler  dissipation  rate  defined  by 


x  =  2D\  V  z 


(4.11) 


The  scalar  dissipation  rate  is  a  very  important  parameter  is  diffusion 
flames.  It  controls  the  mixing  and  determines  the  reaction  rates. 


The  flamelet  equation  for  mass  fraction  is  derived  from  the  scalar  trans¬ 
port  equation  (of  species  k) 


DYk  DYk  d  (  dYk 

Pwyr  +  pUi^ - x—  pD—  I  =  uk 

at  oxj  oxj  \  oxj 


(4.12) 


and  scalar  transport  equation  for  mass  fraction  Z 


dZ  dZ  d  f  „dZ\  n 
p~dt  +  pUid^i  -  d^i  Vd^i  *  “ 


(4.13) 
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Assume  ay  is  a  locally  orthogonal  coordinate  system  attached  to  the  flame 
surface.  x\  is  normal  to  the  flame  surface  and  X2  and  X3  are  parallel  to 
the  surface.  Now  replace  the  coordinate  ay  with  the  mixture  fraction 
Z  and  X2  and  0:3  by  Z2  =  X2,  Z3  =  x3  and  t  —  r.  The  coordinate  Z 
is  normal  to  the  flame  surface  by  the  definition  of  flame  surface.  After 
transformation  of  the  coordinates  ,  Eq.  (4.12)  is  converted  to 


dYk  dYk  dYk\  d(pD)dYk  d(pD)  dYk 
d t  Ui  dZ2  Ui  dZ3 )  dx2  dZ2  dx3  dZ3 


dz\2  d2Yk  dz  d2Yk  dz  d2Yk  d2Yk  &2Yk 
dxj  ~dZ?  +  2  ffx2  8ZSZ2  +  2 dx3  dZdZd  +  dZ\  +  ~dz$ 

(4,14) 


The  flamelct  model  assumes  the  flame  is  thin  in  Z  direction,  the  mass 
fraction  derivative  in  the  flame  surface  normal  direction  is  much  larger 
than  in  the  flame  surface  tangential  direction.  By  an  order- of- magnitude 
analysis  it  can  be  found  that  the  second  derivative  of  Z  and  the  time 
derivative  are  the  dominant  term  on  the  left  side  of  Eq.  (4.14).  Other 
terms  can  be  neglected.  Then  Eq.  (4.10)  for  mass  fraction  of  species 
is  obtained.  The  flamelct  equation  for  temperature  can  be  derived  in 
similar  way. 


To  apply  the  flamelet  model  to  the  simulation  of  turbulent  flames,  it 
assumes  the  “separation  of  scale” :  the  chemical  time  scale  are  small  and 
the  reactions  only  happen  in  a  thin  layer  near  the  flame  surface.  This 
thin  layer  is  assumed  to  have  a  scale  smaller  than  the  turbulence  scale 
and  the  structure  of  the  reaction  zone  remains  laminar  and  unaffected 
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by  the  turbulence. 


One  of  the  widely  used  flamelet  models  is  the  steady  laminar  flamelet 
model  with  a  steady  state  flame  assumption.  The  time  derivative  in 
Eq.  (4.10)  could  be  neglected.  The  governing  equation  of  steady  laminar 
flamelet  model  becomes 


P  Xd2^ 

- — j—  (jJ 

Le  2  OZ 2 


0 


(4.15) 


The  solution  of  the  steady  laminar  flamelet  model  on  Z  space  is  a  function 
of  the  dissipation  rate  and  boundary  conditions.  Thus  the  solutions  can 
be  pre-calculated  and  saved  in  tabular  form.  After  we  solve  Eq.  (4.15) 
and  have  the  solution  of  Z  in  the  flow  held,  the  temperature  and  the  mass 
fraction  of  each  species  can  be  found  by  a  table  look  up.  No  calculation 
is  needed  for  chemical  reactions  during  simulation  since  all  information  is 
precomputed.  This  greatly  reduces  the  time  spend  on  chemical  reactions. 
There  are  also  unsteady  flamelet  model,  such  as  the  “lagrangian  flamelet, 
model”  [55]  and  “Eulerian  flamelet  Model”  [54], 

Larsson’s  simulation  uses  the  “Flamelet /Progress-’ Variable”  model  of 
Pierce  and  Moin  [52],  This  model  uses  the  reaction  progress  variable  C 
as  the  parameter  instead  of  the  scalar  dissipation  rate  y  in  the  flamelet 
model.  The  progress  variable  can  be  defined  as  the  mass  fraction  of  one 
of  the  reactants  or  a  sum  of  several  reactants.  Larsson  [36]  takes  the 
mass  fraction  of  H2O  to  be  the  progress  variable  C .  An  additional  trans¬ 
port  equation  is  solved  for  the  filtered  progress  variable.  Similar  to  the 
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flamelet  model,  the  flamelet  method  uses  precomputed  and  tabulated 
solutions,  parameterised  in  terms  of  the  mass  fraction  Z  and  reaction 
the  progress  progress  variable  C  To  reduce  the  table  dimensions,  the 
operating  pressure  is  assumed  a  priori  as  a  global  constant.  The  mass 
fraction  of  each  species  Yt  and  the  source  term  wh2o  can  be  computed 
as  the  function  of  the  filtered  mixture  fraction  Z,  the  sub  filter  variance 
Z"  Z"  and  the  filtered  C: 


Yi  =  Yi(Z,Z"Z",C)  ,wH20  =  wH2o(Z,Z"Z",C)  ,  (4.16) 

The  filtered  mixture  fraction  Z,  the  sub-filter  variance  Z"  Z"  and  the 
filtered  C  are  calculated  by  the  LES  transport  equations. 

In  Sec.  1.1  we  know  that  in  the  our  model  scramjet  combustor  simulation 
with  high  Reynolds  number,  the  Kolmogorov  scale  (about  10  microns) 
is  smaller  than  the  flame  width,  so  that  the  scale  separation  assumption 
is  not  satisfied.  The  turbulence  eddies  might  penetrate  and  wrinkle  the 
flame  front  and  destroy  the  laminar  flamelet  structure.  Thus  the  locally 
one-dimension  structure  and  might  not  be  a  good  approximation  of  the 
flames  in  the  scramjet  combustion  chamber. 

Instead,  we  use  a  straight-forward  approach  to  simulate  the  turbulent 
combustion  in  our  model  scramjet  combustor  simulation:  we  solve  the 
scalar  transport  equation  for  all  species  for  the  mass  fraction  of  each 
species.  The  chemical  production  rates  are  calculated  by  reaction  rates 
of  each  chemical  reactions  directly.  At  each  time  step,  a  system  of  ODEs 
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based  on  the  reacting  species  are  solved  for  the  change  of  mass  fraction 
of  each  species  from  chemical  reactions: 

=  u.  =  Mk  ^2  VkjTj  (4.17) 

3= 1 

We  refer  to  onr  approach  as  the  “finite  rate  chemistry  approach”  in  other 
parts  of  the  thesis.  The  difference  between  onr  finite  rate  chemistry 
model  and  flam  el  et,  model  is:  we  plan  to  resolve  the  chemical  reaction 
rather  than  model  it. 

Onr  finite  rate  chemistry  model  is  computationally  more  expensive  than 
the  flamelet  model.  Since  we  need  to  solve  systems  of  ODEs  at  each 
time  step  for  every  single  grid  cell.  When  we  use  the  detailed  mecha¬ 
nism,  the  system  of  ODEs  might  include  stiff  equations.  To  reduce  the 
computational  cost  of  this  approach,  a  reduced  chemistry  mechanism  is 
preferred  to  a  detailed  chemistry  mechanism.  We  refer  to  our  approach 
as  the  “finite  rate  chemistry  approach”  in  other  parts  of  the  thesis. 
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4.4  One  Dimensional  Laminar  Flame  Study  with 
Finite  Rate  Chemistry  Model 

In  this  section  we  apply  the  finite  rate  chemistry  model  to  the  simulation 
of  one  dimensional  premixed  and  the  diffusion  flames.  The  purpose  of 
the  chapter  is  to  analyze  the  grid  sensitivity  of  the  finite  rate  chemistry 
model  in  resolving  the  one  dimensional  flames,  to  give  guidance  to  the 
mesh  resolution  required  for  three  dimensional  scramjet  simulation. 

We  considered  two  H2/O2  reaction  mechanisms:  the  detailed  mechanism 
[28]  and  the  reduced  mechanism  [8]  are  considered  here.  We  predict 
resolution  requirements  for  these  two  mechanisms.  We  also  considered 
the  thickened  flame  model. 

4.4.1  Simulation  Configuration 

We  consider  a  long  tube  at  a  constant  pressure.  The  initial  conditions 
of  ID  study  is  in  consistent  with  the  3D  scramjet  simulation  configu¬ 
ration.  For  the  premixed  flame,  H2  and  O2  are  mixed  uniformly  in  a 
stoichiometric  ratio  initially. 

For  the  diffusion  flame,  H2  is  to  the  left,  O2  to  the  right.  O2  has  a 
temperature  of  1500  K,  higher  than  the  H2  flash  point. 

The  hydro  part  of  the  ID  simulation  are  solved  use  the  WENO  scheme 
described  in  Sec.  2.4.  The  molecular  diffusion  coefficients  are  computed 
dynamically  with  the  formula’s  in  Sec.  2.2.3.  Two  H2/O2  reaction  mech- 
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anisms:  the  detailed  mechanism  [28]  and  the  reduced  mechanism  [8]  are 
used  separately  to  calculate  chemistry  source  terms. 


We  try  to  estimate  the  mesh  requirements  for  our  finite  rate  chemistry 
model  to  resolve  the  diffusion  and  premixed  flames.  Here,  a  flame  is  con¬ 
sidered  to  be  resolved  when  certain  measured  parameter  (indicator)  of 
the  flame  is  preserved  and  is  insensitive  to  changes  in  grid  size.  The  lam¬ 
inar  flame  speed  are  considered  to  be  the  indicator  of  premixed  flames. 
The  diffusion  flame  lacks  such  a  simple,  measurable  parameter  as  flame 
speed.  We  use  the  total  heat  release  within  0.12  ms  (one  flow  through 
in  3D  scramjet  simulation)  as  the  indicator. 


4.4.2  Convergence  of  Premixed  Flame 

In  the  numerical  experiment,  initially  we  ignite  the  left  end  of  the  tube 
to  generate  a  flame  transverses  from  left  end  to  the  right  end  in  the  one 
dimensional  tube  with  H2/O2  premixed. 

We  conduct  several  numerical  simulations  of  different  grid  resolution,  the 
grid  size  ranges  from  0.0625  mm  to  1  mm. 

The  convergence  of  the  flame  speed  is  shown  in  Fig.  4.4.  To  achieve 
an  error  of  less  than  10%,  we  need  a  grid  size  of  0.0625  mm  for  the 
detailed  chemistry  mechanism  and  0.125  mm  for  the  reduced  chemistry 
mechanism. 
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Figure  4.4:  Convergence  of  the  flame  speed  of  a  premixed  flame  as  a  function 
of  the  grid  resolution.  The  horizontal  lines  indicate  the  fine  grid  flame  speed 


and  a  variation  of  ±10%  about  this  value.  The  unit  of  the  flame  speed  is  10 


4.4.3  Convergence  of  Diffusion  Flame 

In  the  numerical  experiment,  The  flame  begins  as  the  H2/O2  is  mixed  by 
molecular  diffusion  in  the  middle  of  the  tube.  The  speed  of  the  diffusion 
flame  depends  on  the  diffusion.  The  width  of  a  diffusion  flame  is  generally 
thicker  than  the  width  of  a  premixed  flame.  Thus  the  grid  resolution 
required  to  resolve  a  diffusion  is  less  restrictive  than  the  premixed  flame. 

We  conduct  several  numerical  simulations  of  different  grid  resolution, 
the  grid  size  ranges  from  0.0625  mm  to  1  mm.  The  convergence  of  to¬ 
tal  energy  release  is  shown  in  Fig.  4.5.  To  achieve  an  error  of  less  than 
10%,  we  need  a  grid  size  of  approximately  0.1  mm  for  the  detailed  chem¬ 
istry  mechanism  and  approximately  0.5  mm  for  the  reduced  chemistry 
mechanism. 

The  heat  release  rate  for  a  diffusion  flame  depends  on  the  diffusivity. 
The  convergence  of  heat  release  rate  over  grid  is  also  diffusivity  corre¬ 
lated.  In  2D  and  3D  simulations,  turbulent  mixing  will  speed  up  the 
mixing  process  of  H2/O2  together  with  molecular  diffusion.  We  compute 
the  average  ratio  of  the  turbulent  transport  coefficient  with  the  molec¬ 
ular  transport  coefficient  in  Sec.  4.5.4  and  Sec.  5.2.  Then  we  replace 
the  molecular  transport  coefficient  in  the  ID  diffusion  flames  with  the 
total  transport  coefficient  (sum  of  molecular  and  turbulent  transport  co¬ 
efficients).  The  convergence  of  heat  release  rates  within  1  flow  through 
under  different  grid  resolutions  are  re-examined.  Minimal  mesh  require¬ 
ments  to  achieve  an  error  of  less  than  10%  of  the  heat  release  within  1 
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(a)  The  reduced  chemistry  mechanism.  (b)  The  detailed  chemistry  mechanism. 
Figure  4.5:  Convergence  of  the  heat  release  rate  of  a  diffusion  flame  as  a 
function  of  the  grid  resolution.  Horizontal  lines  indicate  the  energy  release  of 
the  fine  grid  flame  and  a  ±10%  variation  about  this  value. 


flow  through,  for  different  mechanisms  and  flow  regimes  are  summarised 
in  Table  4.2. 


4.4.4  Thickened  Flame  Model 

The  thickened  flame  model  [11]  thicken  the  flame  front  artificially  to 
allow  numerical  solution  of  the  flame  front  on  a  coarser  grid.  The  thick¬ 
ened  flame  model  solves  a  modified  problem,  with  modified  diffusion 
parameters  and  reaction  rate  parameters: 

A  E 

D  ->■  TF  x  D  k  ->•  TF  x  k  u  ->•  TF  x  u  A  ->  — —  E  ->•  — — 

1  p  TF  TF 

(4.18) 
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Here  TF  is  the  flame  thickening  parameter. 


Ying  [72]  has  applied  the  thickened  flame  model  to  the  simulation  of 
one  dimensional  hydrogen-oxygen  premixed  flame.  The  hydrogen-oxygen 
premixed  flame  has  a  very  thin  flame  front,  in  the  order  between  1CT2  cm 
and  10”1  cm  [72],  If  N  cells  (N  >  3)  lines  in  the  flame  width  is  needed 
to  resolve  the  flame,  the  mesh  size  should  be  in  order  of  10-3  cm.  Thus 
resolving  the  premixed  flame  requires  very  fine  grid  resolution.  After 
the  thickened  flame  model  is  applied  the  the  premixed  flame,  a  lower 
resolution  requirement  by  a  factor  of  eight  is  obtained,  in  the  case  of 
a  thickening  factor  of  10  [72],  Thus  the  thickened  flame  model  is  very 
useful  in  the  simulation  of  premixed  flame. 

When  it  comes  to  the  diffusion  flames,  the  thickened  flame  model  delays 
the  ignition  time.  We  conduct  ID  diffusion  flame  simulation  with  dif¬ 
ferent  thickening  factor.  The  heat  release  rate  in  these  simulations  are 
shown  in  Fig  4.6.  We  can  observe  this  ignition  delay  in  Fig  4.6. 

In  the  context  of  our  motivating  simulation  problem,  the  flow  in  the 
scramjet  combustion  chamber  is  supersonic  and  the  residence  time  of 
the  combustible  fuel  in  the  combustion  chamber  is  limited  to  12  ms. 
We  find  that  even  a  modest  level  of  flame  thickening  (TF  =  2),  with 
the  associated  ignition  delay,  reduces  the  overall  heat  release.  Thus  this 
method  is  not  effective.  The  thickened  flame  model  is  excluded  form  the 
model  lists  of  3D  scramjet  simulation. 


(a)  No  thickening 


(b)  Thickening  factor  =  2 


(c)  Thickening  factor  =  4  (d) 


Figure  4.6:  Heat-release  rate  in  diffusion  flame  by  different  thickening  factors. 
The  x  and  y  axis  represent  space  and  time. 
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4.4.5  Conclusion 


We  have  explored  both  the  detailed  and  the  reduced  chemistry  mech¬ 
anisms.  We  find  that  reduced  mechanism  requires  less  stringent  grid 
resolution,  by  a  factor  of  about  2  to  4  in  comparison  with  the  detailed 
chemistry  mechanism.  The  thickened  flame  model  is  excluded  since  it 
lead  to  ignition  delay. 

4.5  Two  dimensional  Turbulent  Flame  Study  with 
Finite  Rate  Chemistry  Model 

The  2D  simulation  is  more  complex  than  the  ID  ones:  it  contains  tur¬ 
bulence,  shock  waves  and  boundary  layers  besides  the  ID  phenomena  of 
molecular  mixing  and  combustion. 

On  one  hand,  turbulence  speeds  up  chemical  reactions.  The  turbulent 
flow  contains  a  wide  range  of  eddies  of  different  length  scales.  These 
eddies  increases  the  mixing  process  the  fuel  and  oxidizer,  allowing  the 
chemical  reaction  to  start.  On  the  other  hand,  combustion  releases  heat 
and  generate  flow  instability  by  buoyancy  and  gas  expansion,  which  in 
turn  enhances  the  turbulence. 

In  this  section,  we  will  study  the  mesh  convergence  of  LES  with  finite 
rate  chemistry  model  in  the  two-dimensional  turbulent  reactive  flow  sim¬ 
ulation,  and  to  understand  turbulent  diffusion,  which  is  also  important 
in  the  three-dimensional  simulation. 
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4.5.1  Simulation  Configuration 


We  choose  our  2D  simulation  as  an  XZ  plane  cut  from  the  model  scramjet 
combustor. 

The  inflow  turbulence  is  generated  by  a  synthetic  turbulent  generator  of 
Tonber  [65].  The  turbulent  intensity  is  2.3%  and  the  length  scale  of  the 
randomised  turbulence  is  about  10  mm.  The  inflow  O2  has  an  average 
density  of  1.241xl0~4  kg/m3,  velocity  of  1890  m/s  and  temperature  of 
1276  K.  The  inflow  H2  has  a  density  of  4.2  x  10-4  kg/rn3,  velocity  of  615 
m/s  and  a  temperature  of  300  K. 

The  grids  for  our  convergence  study  are  defined  by  Table  4.3. 


4.5.2  Energy  Spectrum  Analysis 

The  dynamic  choice  of  parameters  for  the  SGS  construction  compares 
two  resolved  grid  levels,  the  current  grid  and  a  coarser  one,  the  test 
filter.  We  assume  an  asymptotic  (power  law)  behaviour  for  the  dynamic 
coefficient  and  then  the  information  taken  from  these  two  grid  levels  is 
sufficient  to  predict  the  unknown  coefficient  needed  to  model  the  sub¬ 
grid  terms  describing  the  solution  at  a  still  smaller  grid  level.  Thus  the 
fundamental  requirement  for  LES  is  that  the  computational  grid  cutoff 
lies  in  a  scaling  region.  2D  turbulence  has  a  further  complication  with 
two  different  scaling  regions  and  an  inverse  cascade,  those  above  and 
below  the  length  scale  of  the  perturbation  driving  the  turbulence.  From 
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the  above  analysis,  we  take  as  a  convergence  requirement  that  the  grid 
cutoff  must  lie  within  one  of  these  two  scaling  ranges. 

According  to  Kolmogorov’s  hypothesis,  fully  developed  turbulence  has 
three  length  scales  ranges:  the  dissipation  range,  the  inertial  subrange 
and  the  energy-containing  range  [57].  The  inertial  subrange  is  a  univer¬ 
sally  equilibrium  range,  with  an  energy  spectrum  of  E(k)  =  Ce2^3n^5^3 . 
Beyond  the  classical  k~ 5/3  Kolmogorov  scaling,  a  variety  of  other  expo¬ 
nents  have  been  reported  in  a  variety  of  contexts.  Specifically  to  the 
point  here  is  the  2D  nature  of  scaling  exponents  [34,  7],  with  a  classical 
theory  predicting  two  scaling  ranges,  k~ 5/3  for  the  smaller  k  values  and 
k~3  for  larger  k  values.  The  transition  between  the  two  ranges  is  the 
largest  k  value  in  the  turbulence  driver  or  initial  conditions. 

Fig.  4.7  displays  the  temporal  turbulent  kinetic  energy  (TKE)  spectrum 
as  a  function  of  frequency  k  for  several  meshes.  The  data  comes  from  a 
probe  located  in  the  down  stream  of  model  scramjet  combustor  cham¬ 
ber.  We  draw  three  conclusions  from  this  figure.  First,  we  observe  an 
approximate  k~ 5/3  and  k~3  slopes  as  expected.  We  interpret  this  fact 
as  showing  that  the  flow  is  well  resolved.  Second,  the  turbulent  kinetic 
spectra  of  grid  IV  agrees  quite  well  with  the  turbulent  kinetic  energy 
spectra  of  grid  III.  Thirdly,  the  inertial  range,  with  a  slope  of  /c-5//3, 
starts  at  a  frequency  of  2.0  x  105  Hz.  Considering  that  the  average  flow 
speed  is  1500  m/s  at  the  probing  point,  the  grid  levels  I,  III  and  IV 
correspond  to  frequencies  of  1.8  x  106  Hz,  3.6  x  106  Hz  and  7.2  x  106 
Hz  respectively;  the  filter  sizes  of  grid  levels  I,  III  and  IV  correspond  to 
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frequency  (Hz) 


Figure  4.7:  Temporal  turbulent  kinetic  energy  (TKE)  spectra,  with  several 
grid  sizes,  at  a  down  steam  location  in  the  chamber.  The  red  doted  line  shows 
the  At5/3  slope  and  the  black  dashed  line  shows  the  k~3  slope. 
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half  of  these  frequencies.  All  of  these  values  are  strictly  within  the  sec¬ 
ond  scaling  region,  meaning  that  the  SGS  requirement  for  a  grid  cutoff 
within  a  scaling  region  is  satisfied. 


4.5.3  Resolved  Fraction  of  Turbulent  Kinetic  En¬ 
ergy 

Pope  [57]  introduced  a  measurement  of  the  turbulent  resolution  as  the 
comparison  of  the  resolved  turbulent  kinetic  energy  to  the  modeled  ki¬ 
netic  energy.  The  resolved  turbulent  kinetic  energy  is 

kres  =  0.5 ((Hi  -  (ui))(ui  -  (ui)))  .  (4.19) 


The  modeled  turbulent  kinetic  energy  is  estimated  as 


kscs  = 


v2sgs 
Cl  A2 


(4.20) 


where  ^sgs  is  the  turbulent  viscosity  in  the  dynamic  SGS  model  [9].  The 
ratio  M  =  kscs /(kscs  +  kres)  of  modeled  TKE  to  total  TKE  is  used  to 
measure  how  much  turbulent  kinetic  energy  is  modeled  by  the  LES  grid, 
with  smaller  values,  corresponding  to  more  of  the  turbulence  resolved 
rather  than  modeled.  The  spatial  distribution  of  M  values  for  different 
grids  is  displayed  in  Fig.  4.8. 

In  the  interior  of  the  chamber,  we  see  an  increase  in  the  resolved  turbulent 
kinetic  energy.  As  we  refine  the  grid,  we  see  the  higher  M  values  confined 
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grid  I 
grid  II 
grid  III 
grid  IV 


Figure  4.8:  The  ratio  of  modeled  to  total  turbulent  kinetic  energy  in  the 
combustion  chamber 


to  the  boundary  and  inlet  regions.  Thus  the  turbulence  is  well  resolved 
at  grids  II,  III  and  IV  [57]. 


4.5.4  Turbulent  Diffusion 

Turbulent  diffusion  plays  an  important  role  in  modelling  thermal  diffu¬ 
sion  and  the  mixing  of  reactants  at  a  sub-grid  level.  Table  4.4  shows 
the  percentage  of  turbulent  diffusion  coefficient  as  a  fraction  of  the  total 
diffusion  coefficient.  Here  the  turbulent  diffusion  coefficient  are  com¬ 
puted  by  dynamical  Smagorinksy  model  through  Eq  (2.57),  (2.62)  and 
(2.61);  the  molecular  diffusion  coefficient  are  computed  dynamically  as 
in  Sec.  2.2.3.  These  values  are  averaged  over  the  shear  layer  within  one 
scrarnjet  residence  time.  We  observe  that  the  percentage  of  the  turbulent 
diffusion  decreases  with  increasing  grid  resolution. 


In  [45],  a  scaling  law  was  proposed  for  the  turbulent  coefficient  Dcoef  = 
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(a)  The  reduced  chemistry  mechanism,  (b)  The  detailed  chemistry  mechanism. 
Figure  4.9:  Convergence  of  heat  release  rate  of  a  diffusion  flame  as  a  function 
of  the  grid  resolution.  Horizontal  dashed  lines  indicate  the  energy  release  of 
the  fine  grid  flame  and  a  ±10%  variation  about  this  value. 


4 

CAx3.  Columns  4,  5  and  6  of  Table  4.4  show  the  scaled  mean  turbulent 
transport  coefficient.  We  observe  only  a  mild  dependence  of  the  scaled 
coefficient  on  grid  level,  indicating  that  the  scaling  removes  most  of  the 
grid  dependence.  We  re-do  the  ID  convergence  study  of  heat  release 
of  diffusion  flame,  with  new  total  (molecular  plus  turbulent)  diffusivity 
(based  on  the  fraction  of  turbulent  viscosity  out  of  total  diffusivity),  see 
Fig.  4.9.  To  achieve  an  error  of  less  than  10%,  we  need  a  grid  size  of 
approximately  Ax  =  0.52  mm  for  the  reduced  chemistry  and  0.125  mm 
for  the  detailed  chemistry. 
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4.5.5  Reaction  Width  Analysis 


The  width  of  each  reaction  is  shown  in  Table  4.5,  for  a  diffusion  flame 
with  2D  turbulent  mixing.  The  flame  thickness  of  reaction  %  is  defined 
as  [69]: 

f°°  wAx)dx 

Si  =  J-°°  /  u  (4.21) 

maxx\Wi{x)\ 

where  tty  is  the  rate  of  reaction  i.  As  the  turbulent  diffusion  coefficient 
derived  from  the  2D  simulation  is  mesh  dependent,  we  also  specify  the 
mesh  used  to  define  it:  0.125  mm.  We  observe  that  reactions  1,  9,  10 
and  18  have  much  wider  reaction  zones  than  the  remaining  ones. 

For  the  present  analysis,  we  assume  that  at  least  N  (N  >  3)  cells  in  the 
reaction  zone  of  each  reaction  are  needed  for  numerical  resolution.  To 
resolve  the  detailed  chemistry  mechanism  reaction  for  the  2D  diffusion 
flame,  we  need  a  resolution  of  0.22  mm.  To  resolve  the  chemistry  reac¬ 
tion  using  the  reduced  chemistry  mechanism  for  the  2D  diffusion  flame, 
we  need  a  resolution  of  0.34  mm.  Thus  the  reduced  mechanism  has 
a  diminished  resolution  requirement  relative  to  the  detailed  chemistry 
mechanism. 

4.5.6  Grid  Sensitivity  of  the  Chemical  Reactions 

The  reaction  and  turbulence  can  be  observed  from  a  snapshot  of  temper¬ 
ature  and  OH  concentration  in  the  combustion  chamber,  see  Figs.  4.10 
and  4.11.  These  figures  show  the  increase  in  the  resolved  turbulence 


Figure  4.10:  Temperature  snapshot  in  the  combustion  chamber 

grid  I 
grid  II 
grid  III 
grid  IV 

Figure  4.11:  OH  concentration  snapshot  in  the  combustion  chamber 


level  as  the  mesh  is  refined.  The  mesh  convergence  is  measured  by 
two  parameters:  the  total  energy  released  from  the  chemical  reactions  in 
the  combustion  chamber  and  the  percentage  of  H  atoms  that  burn  into 
H2O  [38]  at  the  outlet  of  combustion  chamber.  The  energy  release  and 
percentage  of  H  are  averaged  over  scramjet  residence  time  for  each  grid, 
see  Table  4.6. 

With  an  acceptance  discrepancy  level  set  to  10%,  the  results  indicate 
that  the  combustion  has  reached  grid-convergence  at  grid  III. 
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Table  4.1:  Reaction  Mechanism  of  Hong  [28]. 
Reaction 

Hydrogen-oxygen  chain 

1.  H+02  ^  O+OH 

8.  OH+OH  =  H20+0 

9.  0+H2^  H+OH 

10.  OH+H2^  H+H2Q _ 

Direct  recombination 

18. H2+M^  H+H+M 

h2+h2^  h+h+h2 
h2+o2^  h+h+o2 

7.  H20+M  =  H+OH+M 

H20+H20^  oh+h+h2o 

19.  O+O+M^  02+M 

20.  O+H+M^  OH+M 
Hydroperoxyl  reactions 

2.  H+02(+M)^  H02(+M) 
H+02(+02)^  H02(+02) 
H+02(+H20)^  H02(+H20) 

11.  H+H02^  OH+OH 

12.  H+H02^  H20+0 

13.  H+H02^  H2+02 

14.  0+H02^  0H+02 

15.  H2Q2+H^  HQ2+H2 
Hydrogen  peroxide  reaction 

3.  H202(+M)  ^  20H(+M) 

4.  0H+H202  =  H02+H20 

oh+h2o2  =  ho2+h2o 

6.  2H02  =  02+H202 

16.  H202+H^  H20+0H 

17.  H202+0^  0H+H02 
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Table  4.2:  Summary  of  minimal  mesh  requirements  for  reduced  and  full  chem¬ 
istry,  based  on  ID  flame  analysis,  but  with  laminar  or  2D  or  3D  turbulent 
diffusivity. 


Flow  Regime 

Full  Chemistry 

Reduced  Chemistry 

ID  laminar  flame  analysis 

0.1  mm 

0.5  mm 

ID  flame  analysis;  2D  turbulent  diffusivity 

0.125  mm 

0.52  mm 

ID  flame  analysis;  3D  turbulent  diffusivity 

0.25  mm 

0.625  mm 

Table  4.3:  Grid  resolution 


x  (mm) 

y  (mm) 

Grid  I 

0.5 

0.5 

Grid  II 

0.5 

0.25 

Grid  III 

0.25 

0.25 

Grid  IV 

0.125 

0.125 

Table  4.4:  Turbulent  diffusion  as  a  fraction  of  the  total  diffusion  coefficient. 
The  last  three  column  show  scaled  values  of  these  quantities  to  remove  the 
leading  order  Ax  effect. 


turbulent 

dynamic 

viscosity 

fraction 

turbulent 

thermal 

conductivity 

fraction 

turbulent 
species  (OH) 
diffusivity 
fraction 

scaled 

dynamic 

viscosity 

coefficient 

scaled 

thermal 

conductivity 

coefficient 

scaled 

species  (OH) 
diffusivity 

coefficient 

Grid  I 

0.90 

0.75 

0.73 

0.0000360 

0.00193 

0.149 

Grid  II 

0.28 

0.62 

0.58 

0.0000359 

0.00199 

0.187 

Grid  III 

0.26 

0.40 

0.57 

0.0000326 

0.00203 

0.184 

Grid  IV 

0.16 

0.21 

0.29 

0.0000410 

0.00200 

0.142 
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Table  4.5:  Reaction  widths  for  the  detailed  chemistry  mechanism  in  2D  diffu¬ 
sion  flames, _ 


No 

Reaction 

flame  width  (mm) 

1 

H+02  ^  O+OH 

1.147 

2 

H+02(+M)^  H02(+M) 

1.613 

h+o2(+o2)^  ho2(+o2) 

1.378 

H+02(+H20)^  H02(+H20) 

1.505 

3 

H202(+M)  =  20H(+M) 

0.862 

4 

oh+h2o2  =  ho2+h2o 

1.258 

oh+h2o2  =  ho2+h2o 

1.359 

5 

oh+ho2  =  h2o+o2 

1.467 

6 

2H02  =  02+H202 

0.667 

2H02  ==  02T  H202 

0.890 

7 

H20+M  =  H+OH+M 

1.625 

h2o+h2o^  oh+h+h2o 

1.250 

8 

OH+OH  H20+0 

1.219 

9 

o+h2^  h+oh 

1.100 

o+h2^  h+oh 

1.013 

10 

oh+h2^  h+h2o 

1.096 

11 

H+H02^  OH+OH 

1.440 

12 

h+ho2^  h2o+o 

1.445 

13 

h+ho2^  h2+o2 

1.345 

14 

o+ho2^  oh+o2 

1.417 

15 

h2o2+h^  ho2+h2 

0.959 

16 

h2o2+h^  h2o+oh 

1.054 

17 

h2o2+o^  oh+ho2 

1.313 

18 

h2+m^  h+h+m 

1.325 

h2+h2^  h+h+h2 

1.373 

h2+o2^  h+h+o2 

1.259 

19 

O+O+M^  02+M 

1.319 

20 

O+H+M^  OH+M 

1.240 
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Table  4.6:  Mesh  convergence  of  2D  simulations 


energy  release 

of  reduced 
chem  (J/s) 

energy  release 

of  detailed 
chem  (J/s) 

percentage 
of  H  atoms  in  H2O 

in  reduced  chem 

percentage 
of  H  atoms  in  H2O 

in  detailed  chem 

grid  I 

1.036 

1.106 

21.9  % 

22.3% 

grid  II 

1.392 

1.304 

32.23% 

31.45% 

grid  III 

1.640 

1.541 

36.9% 

34.4% 

grid  IV 

1.775 

1.686 

41.0  % 

40.2% 
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Chapter  5 


Scramjet  3D  Simulation 


In  Chap.  3  we  assess  the  convergence  of  the  averaged  value  of  wall  shear 
stress  and  shear  flux  under  ID  context  and  predict  the  resolution  require¬ 
ment  to  achieve  grid  convergence  within  10%  error  ratio  is  0.25  mm.  In 
Sec.  4.4  and  Sec  4.5  we  made  predictions  regarding  mesh  resolutions 
needed  to  resolve  turbulent  combustion  to  be  0.5  mm.  We  apply  the 
predicted  grid  to  the  simulation  of  3D  scramjet  model  scramjet. 


We  present  a  3D  finite  rate  chemistry  simulation  with  a  mesh  resolution 
of  [Ax,Ay,Az}  =  [0.5  mm,  0.25  mm,  0.25  mm]  for  this  claim.  This  sim¬ 
ulation  compared  with  a  independent  simulation  by  Stanford  PSAAP 
Center  on  the  same  problem,  and  they  are  compared  to  experiment  re¬ 
sults.  The  comparisons  appear  to  be  satisfactory. 
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Figure  5.1:  Sketch  of  the  model  scramjet  combustion  chamber  (Image  source 
[72]). 

5.1  Model  Configuration 

The  scramjet  model  model  combustor  (see  Fig.  5.1)  [21]  studied  here 
is  designed  to  be  geometrically  simple  while  representative  of  realistic 
scramjet  combustion  conditions  in  Mach  5-8  flight  conditions.  The  ex¬ 
periment  setup  is  defined  in  detail  in  [21]. 

The  model  scramjet  combustor  has  an  angled  intake  ramp  followed  by  a 
rectangular  combustor  section  75  mm  wide,  15  mm  high,  and  315.4  mm 
long.  Fuel  is  injected  through  a  single  injector  on  the  center-line  which  is 
D  =  2  mm  in  diameter.  The  distance  of  the  injection  port  downstream 
from  the  flat  plate  leading  edge  is  L  =  70  mm.  Six  high-bandwidth 
pressure  transducers  are  mounted  in  plugs  inserted  into  the  center-line 
of  the  top  wall  of  the  model  combustor  to  allow  pressure  measurements 
[21]- 

The  H2  flows  out  of  the  nozzle  vertically  and  is  bent  downstream  by  the 
cross  flow  02  stream.  The  inflowing  oxygen  has  an  initial  temperature 
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T o2  ~  1200  K,  pressure  Po2  =  0.4  bar  and  cross  flow  velocity  U o2  ~  1800 
m/s.  Together  with  a  kinematic  viscosity  of  v  —  5.36  x  10~4  m2/s, 
the  Reynolds  number  for  the  cross  flow  can  be  computed  as  Reo2  = 
LU/v  =  2.35  x  105,  where  L  is  the  height  of  the  O2  inlet.  The  hydrogen 
is  injected  into  the  combustion  chamber  with  a  pressure  P h2  =  12.5 
bar,  temperature  TH2  =  300K  and  jet  exit  velocity  Uh2  =  1132.5  m/s. 
With  the  kinematic  viscosity  v  =  1.6  x  10~4  m2/s,  the  Reynolds  number 
Re =  DU /v  =  1.4  x  105,  where  D  is  the  diameter  of  the  fuel  jet. 

We  estimate  the  dissipation  rate  e  as  by  e  =  (u'u')3^2/l,  where  l  is  the 
size  of  the  most  energetic  turbulent  eddies  assumed  to  be  a  quarter  of  the 
diameter  of  the  fuel  jet.  The  Kolmogorov  scale  rj  —  (— )  is  estimated 
to  be  approximately  10  microns. 

5.2  3D  Turbulence  and  Turbulent  Diffusion 

3D  turbulence  has  one  scaling  range,  in  contrast  to  the  two  ranges  of 
2D  turbulence.  See  Fig.  5.3,  showing  temporal  TKE  as  in  Fig.  4.7.  Due 
to  the  change  of  the  scaling  exponent  from  -3  to  -5/3,  we  expect  larger 
turbulent  diffusivity  and  a  lower  fraction  of  resolved  turbulence.  In  fact 
the  turbulent  dynamic  viscosity  has  increased  from  2.4  x  10~5  kg/m-s 
(2D)  to  7.3  x  10~5  kg/m-s  (3D),  and  the  fractional  resolution  of  TKE 
plotted  in  a  XZ  cross  section,  Fig.  5.4,  is  to  be  compared  to  Fig.  4.8  and 
still  shows  low  levels  of  modeled  turbulence  within  the  flame  region. 

The  turbulent  viscosity,  turbulent  thermal  conductivity  and  the  turbu- 
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(a)  The  reduced  chemistry  mechanism.  (b)  The  detailed  chemistry  mechanism. 
Figure  5.2:  Convergence  of  the  heat  release  rate  of  a  diffusion  flame  as  a 
function  of  the  grid  resolution.  Horizontal  dashed  lines  indicate  the  energy 
release  of  the  fine  grid  flame  and  a  ±10%  variation  about  this  value. 

lent  species  (OH)  diffusivity  as  a  fraction  of  the  total  viscosity,  thermal 
conductivity  and  species  (OH)  diffusivity  are  approximately  0.65,  0.67 
and  0.84  in  the  3D  simulation. 

We  re-do  the  ID  convergence  study  of  heat  release  of  diffusion  flame, 
with  a  new  total  (molecular  plus  turbulent)  diffusivity  (i.e.,  assuming 
the  same  ratio  between  turbulent  viscosity  and  molecular  viscosity  for 
the  ID  and  3D  models),  see  Fig.  5.2.  To  achieve  an  error  of  less  than  10%, 
we  need  a  grid  size  of  approximately  Ax  =  0.625  mm  for  the  reduced 
chemistry  and  0.25  mm  for  the  detailed  chemistry. 
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Figure  5.3:  Temporal  TKE  spectra  at  a  down  stream  location  in  the  chamber. 
The  doted  line  shows  the  k~ 5//3  slope. 
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Figure  5.4:  The  ratio  of  modeled  to  total  turbulent  kinetic  energy  in  the 


combustion  chamber 
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Turbulence  boundary  layer 


Figure  5.5:  Scramjet  schematic  plot 

5.3  Turbulent  Boundary  Layer  in  Model  Scramjet 
Combustor  Simulation 

The  schematic  plot  of  turbulence  boundary  layer  in  the  model  scramjet 
combustor  is  shown  in  Fig.  5.5.  The  scramjet  simulation  has  solid  wall 
boundary  conditions  at  the  top,  bottom,  front  and  back  walls.  The 
Sec.  5.3  simulation  has  a  grid  resolution  of  (Ax,  Ay,  Az)  =  (0.5  mm, 
0.25  mm,  0.25  mm). 

The  velocity  contour  within  the  combustion  chamber  is  shown  in  Fig.  5.6. 
We  observe  that  the  main  stream  velocity  is  decreasing  in  the  combustion 
chamber  as  the  flow  moves  to  the  outlet.  To  assess  the  boundary  layer 
flow,  the  time  averaged  velocity  normalised  the  main  stream  velocity 
within  the  combustion  chamber  is  shown  in  Fig.  5.7. 

The  boundary  layer  thickness  6  increase  from  inlet  to  the  outlet.  The 
boundary  layer  thickness  £95  at  different  locations  on  each  wall  are  shown 
in  Fig.  5.8.  Here  the  boundary  layer  thickness  in  Fig.  5.8  is  £95,  the  value 
of  wall-normal  distance  at  which  point  the  velocity  reaches  95%  of  the 
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(a)  xz  plane  (y  =  1.6  cm).  Side  wall  is  located  at  y  =  ±3.75  cm 


u:  0  50  100  150  200 


X 

(b)  xy  plane,  (z  =  3.2  cm) 

Figure  5.6:  Velocity  profile  in  the  combustion  chamber 
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x  (cm) 

(a)  xz  plane. (y  =  1.6  cm). 


U/U  inf:  0.4  0.6  0.8  1 


x  (cm) 


(b)  xy  plane,  (z  =  3.2  cm). 

Figure  5.7:  Normalised  velocity  in  the  combustion  chamber.  U-m{  is  the  main 
stream  velocity. 
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main  stream  velocity.  We  plot  £95  instead  of  <599  because  a  1%  velocity 
difference  is  too  small  compared  with  the  velocity  fluctuations  in  the 
interior  part  of  the  turbulent  chamber  to  be  clearly  observable.  Near  the 
inlet,  the  boundary  layer  at  each  wall  has  a  thickness  of  about  1  mm;  At 
most  parts  of  the  top  wall,  boundary  layer  thickness  is  about  3  mm,  At 
most  parts  of  the  bottom  wall,  boundary  layer  thickness  is  about  5  mm. 
At  most  parts  of  the  side  walls,  boundary  layer  thickness  is  larger  than 
7  mm, 


Initially  we  use  the  first  cell  off  the  wall  as  the  matching  point.  Based  on 
analysis  of  Sec  3.5,  an  20%  error  (at  most)  might  occur.  This  potential 
error  explains  the  flow  blockage,  and  the  high  pressure  at  the  outlet.  To 
reduce  this  error,  we  might  need  to  abandon  the  strategy  of  using  the 
first  cell  off  the  wall  as  the  matching  point. 


According  to  the  analysis  of  Sec.  3.5.4,  we  made  a  better  choice  of  the 
matching  point.  The  boundary  thickness  shown  in  Fig.  5.8  means  we 
allow  another  choice  of  the  matching  point.  Near  the  inlet,  we  use  the 
first  point  as  the  matching  point  as  boundary  layer  thickness  5  ~  1  mm. 
At  the  region  x  >  3  mm  we  use  the  interpolated  state  at  y  =  2.5  mm  as 
the  input  to  the  TBL  model.  According  to  error  analysis  of  Sec.  3.5.4, 
the  choice  could  reduce  the  potential  error  in  the  calculation  of  wall  shear 
stress  and  heat  flux  to  under  10%. 


113 


Our  Simulation 

Stanford  PSAAP  Center  Simulaiton 

Turbulent  combustion 

finite  rate  chemistry 

flamelet  /progress- variable 

reduced  mechanism 

detailed  mechanism 

LES 

Subgrid  turbulent  stress 

dynamic  smargorinky  model 

eddy-viscosity  hypothesis  model  of  Vreman 

Subgrid  convective  energy 

gradient  transport  hypotheses 

gradient  transport  hypotheses 

and  scalar  flux 

with  dynamic  Prt  and  Scj 

with  fixed  Prt  =  Sct  =  0.5 

Space  discritization 

uniform  grid 

unstructured  grid,  locally  refinement 

11M  cells 

54  M  cells 

Turbulence  inflow 

Digital  filtering  method 

Turbulent  boundary  layer 

Equilibrium  wall  model 

Table  5.1:  Computational  Set-Ups  of  our  simulation  and  Stanford  PSAAP 


center’s  simulation. 

5.4  Comparison  of  Finite  Rate  Chemistry  and  Flamelet 
Simulations  with  Experiment 

In  this  section,  we  compare  the  finite  rate  chemistry  and  flamelet  simu¬ 
lations  to  each  other  and  to  the  experiments  [21]. 

5.4.1  The  3D  Finite  Rate  Chemistry  and  Flamelet 
Simulation  Models 

Stanford  PSAAP  center  and  we  conduct  large  eddy  simulation  of  model 
scramjet  combustor  with  different  physical  models,  summerized  by  Ta¬ 
ble  5.1.  We  call  our  simulation  “finite  rate  chemistry  simulation”  and 
Stanford  PSAAP  center’s  simulation  “flamelet  simulation” . 

The  finite  rate  chemistry  simulation  uses  a  uniformly  grid  with  11  million 
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cells;  the  mesh  size  is  [Ax,Ay,Az}  =  [0.5  mm,  0.25  mm,  0.25  mm].  The 
flamelet  model  simulation  uses  an  nonuniform  grid,  refinement  near  the 
boundaries  and  the  fuel  inlet,  with  an  overall  mesh  of  54  million  cells. 

The  finite  rate  chemistry  simulation  uses  the  reduced  chemistry  model 
discussed  in  Sec.  4.1,  solving  the  scalar  transport  equation  for  all  species 
and  calculating  the  finite  rate  chemistry  source  terms  directly.  Dynamic 
SGS  models  are  used  to  compute  the  turbulent  diffusion  coefficient  in 
the  LES  transport  equations.  Turbulent  boundary  layers  are  from  the 
method  of  [32], 

The  flamelet  methodology,  explained  in  Sec.  4.3,  is  used  as  the  com¬ 
bustion  model.  The  eddy-viscosity  hypothesis  model  proposed  by  Vre- 
man  [67].  The  sub-grid  heat  flux  and  species  transport  are  modeled 
using  the  gradient  transport  hypotheses  with  fixed  turbulent  Prandtl 
and  Schmidt  numbers.  The  H2  combustion  mechanism  coincides  with 
the  refined  model  referred  in  Sec.  4.1.  Further  details  come  from  the 
related  simulation  [38]. 

5.4.2  Comparison  of  Pressures  on  Upper  Wall 

Fig.  5.9  shows  the  upper  pressure,  comparing  the  finite  rate  chemistry 
simulation,  the  flamelet  model  and  experiment  [21].  Time  averaging 
of  the  point-wise  pressure  removes  turbulent  fluctuations  and  records  a 
mean  pressure  value.  The  smoother  nature  of  the  flamelet  plot  is  due  to 
the  use  of  a  longer  period  for  the  time  averaging. 
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We  divide  the  combustion  chamber  into  the  regions  A  =  [-2,  3]  cm,  B 
=  [3,  7]  cm,  C  =  [7,  11]  cm  and  D  =  [11,  23]  cm.  Region  A  and  D 
show  satisfactory  agreement  among  the  two  simulations  and  the  exper¬ 
iment.  Region  B,  which  is  the  flame  ignition  regime,  shows  satisfactory 
agreement  between  the  finite  rate  chemistry  simulation  and  the  experi¬ 
ment,  while  the  pressure  of  the  flamelet  model  is  too  large.  The  early 
ignition  of  famelet  model  reflects  a  known  weakness  of  steady  flamelet 
models,  which  assume  steady  burning  and  thus  are  unable  to  predict 
ignition.  Consistent  early  ignition  in  “Jet  In  Cross  Flow”  problem  were 
also  observed  by  Chan  [12]  when  she  compared  the  simulation  result  of 
DNS  and  LES  with  flamelet  model.  She  attributed  the  early  ignition  to 
the  omission  of  heat-transfer  into  the  one-dimensional  flame-structure  in 
flamelet  model. 

The  early  ignition  in  the  flamelet-model  can  be  also  be  seen  from  a 
comparison  of  the  water  content  in  region  A  between  the  two  simulations, 
see  Fig.  5.10.  Region  C  is  the  reverse,  with  the  flamelet  simulation 
and  experiment  in  satisfactory  agreement  while  the  finite  rate  chemistry 
simulation  has  too  low  a  pressure. 

5.4.3  Comparison  of  OH  Production 

We  compare  the  OH  concentration  with  the  OH  PLIF  (Planar  Laser- 
Induced  Fluorescence)  plots  from  experiment  at  several  end-view  planes 
from  the  finite  rate  chemistry  simulation,  see  Figs.  5.11  and  5.12.  The 
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Figure  5.9:  Upper  wall  pressure  for  the  finite  rate  chemistry  model,  the  flamelet 
model  and  experiment 
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Figure  5.10:  Instantaneous  snapshot  of  mass  fraction  of  H2O  in  the  scramjet. 
Top  to  bolttom,  frame  (a)  flamelet  (b)  hnite  rate  chemistry,  both  vertical  cross 
section.  Frames  (c)  hamlet,  (d)  hnite  rate  chemistry,  both  a  horizontal  plane. 
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results  appear  to  be  satisfactory,  other  than  the  lower  level  of  OH  in  the 
simulation  along  the  chamber  wall. 

5.4.4  Comparison  of  HoO  Production 

Fig.  5.10  shows  the  water  content  in  the  combustion  chamber  by  simu¬ 
lation  of  finite  rate  chemistry  model  compared  with  the  flamelet  model. 
The  difference  appears  to  be  due  to  two  previously  noted  (see  Sec.  5.4.2) 
features  of  these  simulations.  Namely,  the  delayed  burning  off  the  finite 
rate  chemistry  and  its  slow  burning  after  ignition. 

5.5  Conclusions 

We  see  that  the  3D  turbulent  model  modifies  the  conclusions  of  Secs.  4.4, 
4.5,  by  a  further  increase  of  the  turbulent  diffusion  and  thus  a  further 
increase  in  the  minimal  numerical  mesh  needed  for  convergence. 

Extensive  comparison  of  the  finite  rate  chemistry  and  the  flamelet  model 
simulations  have  been  conducted.  Broadly  speaking,  these  comparisons 
show  consistency  and  substantial  agreement  between  these  two  simula¬ 
tions  as  well  as  consistency  with  experiment.  Most  comparisons  to  ex¬ 
periments  are  qualitative,  as  the  experiments  did  not  calibrate  the  colour 
images  of  concentrations  with  specific  concentration  values.  Figs.  5.11 
and  5.12  show  satisfactory  agreement  of  the  finite  rate  chemistry  sim¬ 
ulation  with  experiment,  with  the  exception  of  a  lower  level  of  chamber 
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x/h  =  11 


x/h  =  13 

Figure  5.11:  Left:  Mass  fraction  of  OH  in  finite  rate  chemistry  simulation 
in  several  end-view  planes  corresponding  to  the  OH  PLIF  data  (right  frame) 
from  experiment 
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(a)  finite  rate  chemistry  simulation 


(b)  flamelet-based  simulation 


(c)  experiment  PLIF  measurement 
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Figure  5.12:  OH  mass  fraction  in  a  vertical  cross  section.,  comparing  finite 
rate  chemistry  simulation,  flamelet  simulation  and  experiment 


(a)  finite  rate  chemistry  simulation  (b)  experiment. 

Figure  5.13:  OH  mass  fraction  in  end  views  of  the  combustion  chamber,  com¬ 
paring  finite  rate  chemistry  simulation  and  experiment 
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wall  burning  in  the  simulation  in  Fig.  5.12.  Fig.  5.10  shows  a  higher 
level  of  burning  for  the  finite  rate  chemistry  simulation,  due  in  part  to 
its  early  ignition,  as  noted  above.  Quantitative  experimental  data  are 
provided  by  the  pressure  measurements  on  the  top  wall.  Agreement  of 
both  simulations  with  this  data  is  satisfactory  with  exceptions  as  noted. 

A  few  consistent  differences  between  the  simulations  have  been  noted. 
The  finite  rate  chemistry  simulation  tends  to  burn  less,  and  starts  burn¬ 
ing  more  slowly.  The  second  of  these  properties,  the  slower  start  to  the 
combustion,  reflects  a  known  weakness  of  steady  flamelet  models  in  not 
capturing  flame  ignition.  Further  details  of  comparison  are  not  mean¬ 
ingful  as  neither  simulation  can  be  regarded  as  fully  converged. 
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Chapter  6 

Conclusions  and  Future  Work 

6.1  Conclusion 

In  this  thesis,  we  study  the  physical  processes,  including  the  turbu¬ 
lent  boundary  layer,  turbulent  mixing  and  turbulent  combustion,  in  the 
model  scramjet  combustor  using  Large  Eddy  Simulation. 

In  Chap.  3  we  explore  the  turbulent  boundary  layer  model  and  assess 
mesh  convergence  requirements  for  the  averaged  value  of  wall  shear  stress 
and  heat  flux  through  suits  of  two-dimensional  simulations. 

In  Chap.  4,  we  examine  the  mesh  requirement  for  the  finite  rate  chem¬ 
istry  model  to  resolve  the  turbulent  combustion  through  suits  of  one- 
dimensional  and  two-dimensional  simulations.  There  are  three  predic¬ 
tion  approaches: 

(a)  Convergence  study  of  one-dimensional  laminar  diffusion  flame.  We 
first  conduct  the  simulation  with  laminar  diffusion  coefficients.  Later 
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Table  6.1:  Summary  of  minimal  mesh  requirements  for  reduced  and  full  chem¬ 
istry. 


Flow  Regime 

Full  Chemistry  (mm) 

Reduced  Chemistry  (mm) 

ID  laminar  flame  analysis 

0.1 

0.5 

ID  flame  analysis  with  2D  turbulent  diffusivity 

0.125 

0.52 

ID  flame  analysis  with  3D  turbulent  diffusivity 

0.25 

0.625 

2D  flame  width  analysis 

0.22 

0.34 

we  add  the  turbulent  diffusion  coefficient  from  2D  and  3D  simula¬ 
tion  onto  laminar  diffusion  coefficient  and  get  a  less  stringent  crite¬ 
ria. 

(b)  Two-dimensional  flame  width  analysis. 

The  results  are  summerized  in  Tab.  6.1.  Top  four  row  of  Tab.  6.1  predict 
the  mesh  needed  for  the  finite  rate  chemistry  model  to  resolve  the  chemi¬ 
cal  reactions.  They  show  some  level  of  consistency  (although  not  perfect 
consistency).  We  also  find  that  reduced  chemistry,  which  removes  the 
fast  and  transient  reactions,  has  a  less  stringent  mesh  requirement  than 
the  full  chemistry. 

In  Chap.  5  we  applied  the  predictions  made  in  Chap.  3  and  Chap.  4  to  the 
3D  simulation  of  model  scramjet  combustor.  TKE  analysis  and  turbulent 
resolution  analysis  in  3D  scramjet  simulation  shows  the  turbulence  in 
them  combustion  chamber  is  well  resolved. 

The  simulation  result  are  compared  with  the  experiment  and  Stanford’s 
simulation.  Qualitative  comparison  of  OH  mass  fractions  show  satisfac¬ 
tory  agreement  of  the  finite  rate  chemistry  simulation  with  experiment, 
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with  exceptions  of  a  lower  level  of  chamber  wall  burning  and  a  higher 
level  of  burning  near  outlet.  Quantitative  measurements  of  top  wall  pres¬ 
sure  shows  satisfactory  agreement  of  both  simulations  with  experiment 
except  for:  flamelet  simulation  has  a  higher  level  of  burning  due  to  early 
ignition;  finite  rate  chemistry  model  has  lower  level  of  burning  due  to  ig¬ 
nition  delay.  The  satisfactory  agreement  of  simulation  with  the  Gamba’s 
experiment  data  verifies  predictions  made  in  Chap.  3  and  Chap.  4. 

6.2  Recommendation  for  Future  Work 

Although  the  current  scramjet  simulation  result  presented  in  this  thesis 
have  demonstrated  good  agreement  with  the  experiment,  it  could  be 
further  developed  in  a  number  of  ways. 

First  recommendation  is  the  implementation  of  local  mesh  refinement 
functionality  into  current  FronTier  package.  The  flame  thickness  list 
in  Tab.  4.5  is  the  averaged  flame  thickness  in  the  combustion  chamber. 
Actually  the  flames  near  the  H2  inflow  nozzle  has  a  smaller  width  than 
the  averaged  value.  So  it  is  better  to  have  local  mesh  refinement  near 
H2  inflow  nozzle.  We  can  also  have  less  error  in  the  calculation  of  the 
wall  shear  stress  and  heat  flux  with  mesh  refinement  near  wall. 

Another  recommendation  is  the  time  efficiency  optimization  of  FronTier 
code  especially  in 

(a)  Parallel  Computing. 

The  parallel  computing  of  scramjet  project  is  realized  by  MP1.  In 
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the  implemenation  of  parallel  computing,  the  scramjet  simulation 
domain  are  split  into  several  uniform  rectangles  (subdomains)  and 
each  processor  is  in  charge  of  one  subdomain.  This  kind  of  splitting 
is  not  time  efficient  because  the  computational  burden  is  not  evenly 
distributed  on  these  subdomains.  Some  subdomains  contains  lots 
of  “obstacle  states”  and  are  fast  than  others.  Thus  lots  of  compu¬ 
tational  resources  are  wasted  when  the  fast  subdomains  stops  and 
waits  for  the  slow  subdomains.  Thus  a  better  subdomain  splitting 
algorithm  are  needed. 

(b)  Computation  of  EOS  and  transport  coefficients. 

In  scramjet  simulation,  EOS  and  transport  coefficient  functions  are 
called  quite  frequently  by  hyperbolic,  parabolic,  boundary  layer  and 
chemistry  modules.  Calculation  of  the  EOS  parameters  and  the 
transport  coefficients  for  gas  mixtures  require  lots  of  computational 
resources.  Optimization  of  the  EOS  and  transport  coefficients  code 
will  greatly  reduce  computational  cost. 
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Appendix  A 

Mathematical  Derivation  of  the  Thin 
Boundary  Layer  Model  for  Compressible  Flow 


In  this  section  we  shows  the  derivation  the  Thin  Boundary  Layer  Model  of 
[4]  from  the  Navier-Stokes  equation.  The  compressible  Navier-Stokes  equation 
expresses  conservation  of  mass,  momentum,  and  energy.  The  law  of  conserva¬ 
tion  of  mass  is 

tr  +  -^r(W)=0-  A-1) 


m  +  d^[pu0=0' 

Conservation  of  momentum  is  given  by  the  equation 


d  (  \  d  (  \  dp  dTij 

at  +  fe  =  ' 


(A. 2) 


The  equation  for  conservation  of  energy  is: 


dE  d\(E  +  p)  uj]  d 

dt  dx-j  dx.j 


(A.3) 


For  a  Newtonian  fluid,  assuming  Stokes  Law  for  mono-atomic  gases,  the  vis¬ 
cous  stress  is  given  by: 


_ r)  o*  _  \fdiii  duj\  2  duk' 

Tij  llSv  11  Idx^dxJ  3 5ijdxk  • 


(A.4) 


The  energy  equation  could  also  be  written  in  terms  of  physical  variable  hb: 


dhb  dhb 
dt  Ul  dxj 


dp  d  dq3 

dt  +  d.r,  dxj 


(A.5) 
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In  this  equation,  hb  is  defined  as: 


hb  =  h+  \ul  ,  (A. 6) 

and  the  specific  enthalpy  is  given  by 

h  =  e  +  -  (A. 7) 

P 

where 

e  =  cvT  (A. 8) 


is  the  specific  internal  energy.  Some  authors  define  hb  to  be  the  total  enthalpy 
and  h  as  the  static  enthalpy.  We  assume  the  perfect  gas  law  is  applicable  here, 
which  assumes  that 


Mr  =  <7-1KT, 

(A.9) 

h  =  cvT  +  (7  -  1)  cvT  =  7 cvT  =  cpT  . 

(A. 10) 

Thus 

hb  —  cpT  +  —  uk  . 

(A. 11) 

The  heat  flux, 

q  is  given  by  Fourier’s  law 

q  =  —kVT . 

(A. 12) 

Within  the  turbulent  boundary  layer,  we  replace  the  instantaneous  quantities 
in  the  Eqs.  (A.l  -  A. 3).  by  the  sum  of  their  mean  and  fluctuating  parts.  We 
introduce  the  Reynolds  decomposition  for  p  and  p  and  the  Favre  decomposition 
for  T,  Ui,  hb  and  h.  By  taking  an  average  of  the  governing  equations  we  have: 


dp  d 
dt  dxj 


0, 


d  1 

\  d  , 

dt\ 

pUj 

\ 

/  dxj ' 

d 

- 

d  r 

di 

phb 

+  dxj  [ 

pUiUj 


dp  +  dr  jj 
dxi  dxj 


pUjhb  +  pu-h'l  -  TijUi  +  qj 


ay  H'O 

_  dP 
~  ~dt  ’ 


(A.13) 


(A.  14) 
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where  the  density  averaged  hb  is  give  by 


~  ~  UkUk  P  u , 

hb  =  cvT  +  +  -  +  -f  . 

a  p  z 


(A.15) 


Next  we  substitute  the  unknown  terms  into  the  known  terms.  For  example  we 
have 

Tjt  =  TZ  +  Tji,  (A.16) 

'  ^  (A. 17) 

(A. 18) 


//  7  n  /  rr,  -1-  o\  'f  rj-\  .  //  //  .  -1-  n  A  n 

pUjh  =  puj  (cpT  +  -uf)  ==  CppuJP  +  uipuiuj  +  -pu  %u3  , 
_  ,  dT  r  8 f  ,  9T" 

_  _A:a^” _  -A;a^“  “  a^“  ’ 


and 


UjT ij  XliTij  T  W j  Tij  T  Vj^Tjj  . 

Then  we  rewrite  equations  Eqs.  (A.l  -  A. 3)  as  : 

dp  d  f 


a  /  _\  d 

m  + 


ax,- 


+  Pui  uj  ~  Pj  -  Tij  }  =  dx,i 


dp 


d  d  _  _ 

at  +  ( ' ?Ujhb  +  Uipu'iu'i +  ^ Pu'juiui 


a 


ax,- 


-cppu,T  -  k 


dT 

dxi 


dT" 


UiTij  k 

ax,- 


r*i  Ta  ]  m 


(A. 19) 

(A. 20) 
(A. 21) 


=  •  (A-22) 


These  are  open  equations  in  that  they  contain  new  undefined  terms.  Clo¬ 
sure  requires  same  assumptions.  The  Reynolds  stress  r^rb  =  pu”u"-  can  be 
modelled  using  an  eddy-viscosity  assumption: 


_ turb 


nr" =  -puiuj  ~ 2 TtS*j  -  -pkdij , 


(A. 23) 


where  pt  is  the  turbulent  viscosity.  The  term  cppUjT  corresponding  to  turbu¬ 
lent  transport  of  heat,  can  be  modelling  using  a  gradient  approximation  for 
the  turbulent-flux:  _ 

Pt  dT 

<1.1 .  "  '  “ 


Jurb  =  cpPu”  T  =  cf 


Prt  dxj 


(A. 24) 
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Here  Prt  is  the  turbulent  Prandtl  number,  which  is  taken  to  be  a  constant 
with  value  0.9.  The  term  rrj  is  neglected  if  |  Ay  |>>|  |  which  is  true  for 

nearly  all  flows.  The  terms  cpk^~,  pUj  u'  and  rTaU-  are  neglected  for  similar 
reasons.  Thus  the  filtered  averaged  equations  for  compressible  flow  become: 


dp  d 

m  +  W,  K1  =  0 ' 

pui  d  r__  __  f  dui  duj  2  duk 

at  +  &/  b“’  - ("  +  "l)  a&A 


(A. 25) 


+  iH'  (A-26) 


dt  dxj 


phb  d  _  ~  ~  (  .  v  ( dui  duj  2duk  \  /  p  \  dT  _  dp 

m  Pujhb  Ui  (/i  +  /it)  o  ~  o  cm  +  (  cp  +  ^ 


<9ay  <9ay  3  dxk 


(A. 27) 

For  the  boundary  layer  flow,  we  use  an  order-of-magnitude  analysis  to  derive 
a  simplified  system  of  equations.  If  y  is  the  wall  normal  direction,  we  assume 
that 

„  _  __  d  d  d 

V<£u,w,  — ,  —  <  — 
ox  oz  ay 

Omitting  lower  order  terms,  we  have  the  new  boundary  layer  equations  (for 
three  dimensional  flow): 


T  (pB)  +  T  (p5)  +  =  o  , 

pu  _~du  _~dw  _~du  dp  d  du 

m  +  ^  +  &  =  %  +  ^  dy 


(A. 28) 


(A. 29) 


phb  _~dhb  _~dhb  _~dhb  d\~,  \du  pt  dT  \  1  dp 

— +pu  — +  pu  — Tprc— =  -^^  +  ^-+^  +  0,  — 

(A. 30) 

Our  model  is  a  equilibrium  stress  model  without  pressure  gradient  so  that  the 
left  hand  side  of  the  momentum  equation  is  zero.  According  to  Bernoulli’s 

u2 

law,  hb  =  h  +  -y  along  a  stream  line  is  constant,  so  the  left  hand  side  of  the 
energy  equation  is  zero.  So  our  boundary  layer  equation  is  further  simplified 
to  become: 

d  T  du  1 

yr  (h  +  Pt)  yr  =0,  (A.31) 


/it  dT 


131 


d_ 

dy 


dT 

dy 


+  (/^  +  A*t)  u 


du 

dy 


0  . 


(A. 32) 
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